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ABSTRACT 


Advances in technology have given way to concepts in warfare that were once 
constrained to the world of science fiction. In an effort to stay ahead of any potential 
adversary’s weapons development, we must look down the path of countermeasures to 
high-energy electromagnetic weapons. 

In the attempt to engineer a material that can reduce transmitted beam intensity by 
the greatest factor, we look to liquid crystals. They have great potential to provide a 
starting point to engineer a material in order to show increased protection of DoD assets 
from high-energy beam weapons. 

We first develop one-dimensional finite-difference time-domain codes to solve 
Maxwell’s equations in order to model the electromagnetic wave propagation in a liquid 
crystal layer. After validating numerical results with analytical results for matched 
anchoring, we investigate the heterogeneous liquid crystal structures with mismatched 
anchoring conditions and determine the best anchoring conditions to minimize 
transmitted beam intensity. The main result of the simulation is that for a known incident 
wave the maximum reduction of the transmitted intensity is achieved with matched 
anchoring conditions. However, for mixed anchoring conditions, there is evidence that 
the mixed structure can reduce the intensity for a wider range of waves. 
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I. ELECTROMAGNETIC WAVE PROPAGATION 
IN LIQUID CRYSTALS 

A. INTRODUCTION 

War is a “come as you are” business; on the onset of conflict, we can only 
bring what we have and only until we realize the capabilities of our enemy, can we 
then decide to try and develop means by which to counter their offensive. Therefore, we 
must stay ahead of any potential adversary’s weapons development and anticipate 
what could be brought to the table, should that conflict occur. The world of high-energy 
electromagnetic weapons is no different. 

It is through this need to stay ahead that we must look at the possibilities of 
utilizing known material or engineering a material with the capabilities to minimize the 
impact of these types of weapons. For this, we look to liquid crystals (LCs). 

Through the continued study of propagation of electromagnetic energy in liquid 
crystals (Hwang & Rey, 2005; Hwang, Han, & Rey, 2007; Kriezis & Elston, 1999; 
Kriezis & Elston, 2000; Ziogos & Kriezis, 2008), it becomes increasingly important to 
branch off of the typical models that utilize a consistency in alignment of anisotropic axes 
with that of the coordinate axes which has been studied and looked at by (Choate & 
Zhou, 2012.) In order to push passed the constraints of that study, we must first validate 
the model used by Choate and Zhou and then use the validated model to investigate the 
differences between matched and mixed anchoring conditions of the liquid crystal layer. 
This shifts our focus from the matched anchoring conditions on the top and bottom of the 
homogeneous layer to mismatched or mixed anchoring conditions of an anisotropic layer. 

When considering liquid crystals, we see that there are many advantages to 
utilizing them in an effort to minimize or maximize the intensity of a transmitted beam. 
For the purposes of this paper, we will take a generic look at liquid crystals but with the 
specific focusing on the way that the molecules within the liquid crystal medium align. 
Additionally, even though there have been studies about maximizing the transmitted 
intensity, we will only be looking at efforts to minimize that transmitted intensity. 
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In Figure 1, we show a cartoon of the problem setup. We consider a 
monochromatic electromagnetic wave passing through a liquid crystal layer. The liquid 
crystal layer is sandwiched between two supporting glass plates. It drives the question as 
to how adjusting the internal structure of a liquid crystal layer impacts the difference 
between the intensity of an incident beam as compared to the intensity of the transmitted 
beam through the layer. Our ability to adjust the internal structure of the material lies 
within the anchoring conditions at the top and bottom plates located at each end of the 
layer. In this thesis, we assume that the orientation of the liquid crystal molecules is not 
affected by the external electromagnetic field. 


Glass plate 


/// /// '/"A' 


/// 

a 


LicuieCryv.^P la 


Gneesisoe uutned » be fried tea 
ttr. effected fey ae 5AS 



Anchoring conditions induce a 
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propagation"’ 


Glass plate 


Incident Beam 


Figure 1. A cartoon description of the problem setup. 

B. A BRIEF OVERVIEW OF LIQUID CRYSTALS 

The concept of liquid crystals was realized in 1888, by the Austrian chemist 
Friedrich Reinitzer. While engaged in studies with regard to cholesterol, he noted that 
the material cholesteryl benzoate had what appeared to be two melting points. The first 
melting point at 145.5 degrees Celsius resulted in a cloudy liquid that remained cloudy 
until heated to the next melting point of 178.5 degrees Celsius. At this higher 
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temperature, the material became a transparent liquid. Further investigation revealed no 
impurity based cause for this phenomenon. 

Through efforts of Reinitzer and crystal optic expert and German physicist, Otto 
Lehmann, it was determined that the cloudy state was unique and presented a material 
that was between a liquid and solid with shared properties of both. This opened the door 
to the world of liquid crystals, a material with anisotropic properties that varied with 
orientation of molecules (Liquid Crystals, 2013.) Solids have a required order of both the 
orientation and position of the individual molecules of the substance. Liquids, on the 
other hand, are not bound by this same principle and can change the orientation and 
position of the molecules. The discovery by Reinitzer and Lehmann revealed that the 
material, termed liquid crystals, had an orientation order but no position order. This 
explains why they chose to term the material as a liquid crystal. 

It is through these varying anisotropic properties that we choose to utilize the 
concept of a liquid crystal as the medium by which to affect transmitted beam intensity. 
One method used to effect the orientation of the molecules within the liquid crystal is to 
establish the anchoring conditions on a surface by either a mechanical rubbing or 
grooving of the surface. This causes the molecules to align with the mechanical aspect of 
the grooves. Another method is to apply an electric charge to the layer, thereby effecting 
orientation. Once orientation is established by anchoring conditions, the other molecules 
in the immediate vicinity of those anchored molecules will align their orientation, with 
respect to y/ and 0 , similar to that of the anchored molecule as illustrated in Figure 2 and 
Figure 3. With that in mind, we must establish that the orientation is known and 
predetermined prior to the incident beam interaction and that the electromagnetic field of 
the beam will have no impact on preexisting molecule orientation. 

We use the Leslie-Ericksen (LE) model theory to help us describe the behavior of 
the liquid crystal which has the characteristics of a rigid-rod system. For a more general 
description of the mathematical modeling aspects of liquid crystals, see (Zhou, Forest & 
Wang, 2010.) As seen in Figure 2, the molecules have a high aspect ratio, which in turn 
allows for anisotropic characteristics. The molecule’s centers of mass are random but 
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there is an orientation order coinciding with the individual molecules. Based on LE 
theory, we define the unit vector m as the axis of symmetry of a spheroid molecule 
where m is not random. We further define the unit vector n as the major director that is 
the preferred direction of m at the given (x,r) . The major director is the local average 
orientation of molecules. 



Figure 2. Nematic phase characteristics of LCs. 


The nematic phase nature of liquid crystals means that the molecules of the 
material possess orientation direction for the given property conditions but requires no 
positional direction and can therefore be considered random. 

The orientation direction of the material is described by a major director axis. 
Through the assumptions that no reorientation of molecules occurs by the 
electromagnetic field of the incident beam or by the flowing of material within the liquid 
crystal, we can establish and utilize the major director, nas determined by the steady 
state equation 

(I-nn)-V 2 n = 0 (1.1) 


4 


where nn is the dyadic product, V = 


d d d 
dx dy dz 


"I t 


is the gradient operator and 


V 2 = V • V = —--t-— H-- is the Laplacian operator. That is, V 2 n is in the null space 

dx dy~ dz 

of I - nn , which is An , where A is an arbitrary scalar. Therefore, we need to solve 

V 2 n = An . (1.2) 

We are left with three unknowns and an arbitrary A that we can choose in order to solve. 

Since n is a unit vector, we can rewrite n in terms of the angles 0 and (// as 

cos 6 cos y/ 


n 


cos 6 sin y/ 
sin# 


(1.3) 


Plugging this form into Equation(l.l), we have a set of equations for both 0 and (//, 

V 2 y/ - 2 tan (N 0 ■ V y/ = 0 
V 2 # + sin#cos#V y/ ■ Vy/ = 0 . 

V G • V 6 + cos 2 (N y/ -V y/ = —A 


We only need to solve the top two equations, with the understanding that A is arbitrary 
and its actual value is irrelevant, thereby reducing the system of equations to 


V 2 ^-2 tan #7#- V y/ = 0 
V 2 # + sin#cos (Nyr -Wy/ = 0 


(1.4) 


The solutions to these equations are determined by the boundary conditions. In 
our specific case, these boundary conditions are established by the two glass plates which 
support the liquid crystal layer at z = 0 and z = h. These plates will be referred to as the 
bottom and top plates, respectively. Utilizing anchoring conditions that are uniform in the 
x - and y -directions, the system of Equations (1.4), can be reduced to a system of second 
order ordinary differential equations in the z -direction. In order to solve, this requires 
four boundary conditions defined as 
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o hol = e{z = 0) 
v b o,=v(z = Q) 

0, op = 0{z = h) ■ 

V 7 to P = viz = h) 

These four boundary conditions are our control parameters for determining the 
orientation of molecules within the layer and the resultant propagation of light through 
the layer. The goal of this paper is to determine the effect of varying these parameters on 
the transmitted intensity of a given incident beam. As seen in Figure 3. , we define 6 as 
the angle between the xy - plane and the z — axis. The value of 0 is limited to values of 
0 to k 1 2 where 0 = 0 is at the xy - plane. Additionally, we define <// as the angle that 
lies within the xy - plane. The value of y/ is limited to values of 0 to n / 2 where 
y/ -0 is in line with the x- axis. 



y 


Figure 3. Description of 6 and <// . 

It is through the values of 9 and <// at the boundaries that we can effectively tune 
the layer to achieve maximum or minimum intensity of the transmitted beam. We will 
show later through model runs of how to best determine 6 and t// for our given 
assumptions. 
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C. MAXWELL’S EQUATIONS 

This study will begin by first reviewing the fundamental basics of previous 
studies of electromagnetic plane waves and their propagation through liquid crystal 
medium of constrained orientation. The orientation structure is determined by anchoring 
conditions on the top and bottom plates. Our goal is to find the dependence of the 
transmitted intensity on the structure, controlled by the anchoring conditions. Following 
that, we will validate the model used in previous work (Choate & Zhou, 2012) through 
the comparison of incident beam intensity to the transmitted beam intensity across the 
given medium. Finally, we will mn a modified model that will vary the angles of 
orientation within the given medium to determine overall impact of the transmitted 
intensity of the given energy beam. 

As outlined in work done by Choate and Zhou (Choate & Zhou, 2012), we will 
maintain consistency in the coordinate system and geometry of the overall system. This 
choice will allow for a more convenient comparison of data provided by their past studies 
and it will also allow the validation of the model used in their studies. In order to derive 
the equations that we will be working with, we must first establish the layer type that we 
will be working with. 

As shown in Figure 4, the regions z < 0 and z > h are vacuum regions with the 
given permittivity and permeability, s o and //„. The region between z = 0 and z = h , is 
the medium of focus and is modeled by a liquid crystal polymer of thickness h . The 
orientations in the x - and y -directions are assumed to be homogeneous but are not 
required to be the same value for v as for y . The orientation is given by the major 
director, represented by the unit vector n, that can vary in space. It is determined by 
Equation (1.1) coupled with the anchoring conditions at the plates. Since the anchoring 
conditions are uniform in x- and y-directions, we have a problem. We want to examine 
the effect of using different anchoring conditions on the top and bottom on the intensity 
of the transmitted beam. Later, we will solve for the intensity of the transmitted beam. 
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z=0 z=h 


Figure 4. Liquid crystal layer orientation and geometry. 
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The fundamental starting point for any electrodynamics system is to look at 
Maxwell’s equations 


— = VxH (1.5) 

dt 


V D = 0 


( 1 . 6 ) 


SB 

dt 


-VxE 


(1.7) 


V -B = 0 (1.8) 

where E is the electric field, D is the electric displacement field , H is the magnetic field, 
and B is the magnetic flux density. 


We assume that the liquid crystal is nonmagnetic, such that B = // 0 H , where // 0 is 

the permeability in free space. The anisotropy of the liquid crystal appears in the 
constitutive equation for the electric field, 

D = £-E (1.9) 

for the permittivity tensor, 
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£ = fJ+Affim (1.10) 

where a s = e n -e L . Solving the system of Equations (1.4) for (// and 0 and substituting 

them into Equation (1.3) gives us the n that can be substituted into Equation (1.10) for the 
permittivity tensor. The extraordinary permittivity, < 5 j |; is the permittivity in the direction 

of n . The ordinary permittivity, £ , is the permittivity in directions orthogonal to n . We 
will assume that > £ due to the rod shape of the liquid crystal molecules. 


We introduce the rescaling, 

d 1 




D, E = ,P-E, l = ± 
Mo 


£ o Mo 

where e 0 is the vacuum permittivity. While the relative permittivity £ is dimensionless, 

both D and E are not dimensionless. They have equivalent dimensions to H. Also, the 

1 


free space quantities £ 0 and ju 0 are defined such that 


light in a vacuum. This rescaling transforms the system (1.5)(1.7)(1.9) into 

1 


= c 0 , where c 0 is the speed of 


dD 


^ \J £ oMo 


VxH 


D = £ E 
5H 1 


dt 


yJ £ oMo 


VxE 


( 1 . 11 ) 

( 1 . 12 ) 

(1.13) 


Since the anchoring conditions in the x - and y -directions are homogeneous, the 
partial derivatives are zero in those directions. Utilizing the condition that all partial 
derivatives with respect to x and y are zero, the components of Equations (1.11) 
through (1.13) are: 


dD x 

dt 

ii 

i 

O 

(1.14) 

8D y 

dt 

O 

II 

(1.15) 
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_ QJ 

II 

o 

(1.16) 

D = g E, 

(1.17) 

o 

II 

CO 

(1.18) 

dH y dE x 

dt C ° dz ’ 

(1.19) 

8H < =0. 
dt 

(1.20) 


The form and nature of Equation (1.16) and Equation (1.20) tells us that the 
components D_ and H_ are constants with respect to time. Knowing that the H and /) 

components of the incident beam are zero, we set them as H_ — 0 and D_ = 0 through 
the remainder of the layer. However, E. is nonzero inside the layer. As a consequence of 
D_ = s xz E x + s E + s,,E. = 0 combined with Equation (1.17), we have 

z ~ v xz X yz yS 

G zz 

Once the electric and magnetic fields are known, we can compute the intensity as 
follows. First the Poynting vector is defined as the cross product of the electric field and 
the magnetics field 



The intensity is defined as the time average over one period of the magnitude of the 
Poynting vector or 



For convenience, we drop the tildes for the remainder of the paper. 
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D. NUMERICAL ALGORITHM 

We will solve the system. Equations (1.14) through (1.20) numerically utilizing 
the one dimensional finite-difference time-domain (FDTD) method (Taflove & Hagness, 
2005.) This method uses the Yee cell discretization where the electric field and the 
magnetic field are stored on staggered grids. This allows for a second order scheme in 
both space and time while only requiring storage of one field value per grid cell. The 
problem that we are trying to solve assumes an infinite domain in the z -direction, 
therefore we have to truncate the computational domain with perfectly matched layers 
(PMLs) to minimize the artificial reflections created by the truncation of the domain with 
perfect electric conductor (PEC) boundaries. 

The incident beam is a plane wave where the wavelength, polarization, incident 
angle, electric field E mc , and magnetic field H /)ic are known. These fields represent the 

beam that would pass through the computational domain if there were no liquid crystal 
layer and are assumed to be known throughout the computational domain. The incident 
angle is supposed to be perpendicular to the bottom plate and in line with the z — axis. In 
order to introduce the incident beam to the model, we split the computational domain into 
a total field (TF) region and a scattered field (SF) region, as shown in Figure 5. We 
define the scattered fields E scaf and H scaf as the differences between the unknown fields 
and the known incident fields so that the total unknown fields are 

E.., = E +E. , H,,, = H +H . 

In the TF region, we store and update the total fields, but in the SF regions, we 
only store and update the scattered fields. This allows us to impose the incident beam at 
the z start in the SF region without the nonphysical reflections that would be caused by the 

reflected beam returning to z start . This convention does introduce the difficulty of 

calculating the electric and magnetic fields as they pass from one region to another, 
which we will discuss later. 
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Scattered Scattered 



- Insertion of the incident beam at Z start 

Figure 5. Computational domain. 

Our computational domain consists of the liquid crystal layer between the points 
z glassbot and z Rlmstop ( bein g z gIassbot + h) which is symmetrically surrounded on both ends 

with a region of vacuum between the points z glasstop and z PMLtop and a PML region 

between the points z PMLtop and z PEChot , which marks the end of the truncated domain. 

There exist similar layers beyond the bottom end of the liquid crystal layer as shown in 
Figure 5. Additionally, between z glasstop and z PMLtop , we have a transition between the 

total field to the scattered field at the point z scattop and similarly at z scathot for the bottom 

portion. The incident beam is introduced in the bottom scattered field at z stan located 

between z PMLhot and z scathot . 

In order to discretize space and time, we define the grid spacing by the width of 

h 

the liquid crystal layer, h , between the bottom and top plates, where az = —, and N is 

the number of cells within the layer. We define the grid points z. k = k&z for 
kpEctot - k - kpEcmp where k- 0 corresponds to z glassbot . Similarly, for time, t n = ha t for 
n = 0,1,2... where a t is determined by the FDTD method. The value of t fmal is sufficient 
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such that the transient wave can propagate completely across the liquid crystal layer with 
two additional wave periods to allow for averaging to compute the intensity of the 
transmitted beam. 


The electric field components are computed to coincide with the grid points. 



(1.21) 

D;, k =D y (z k ,t n ) 

(1-22) 

Kk= E MM 

(1-23) 

E n y , k = E y (z k ,t n ). 

(1-24) 


In order to maintain the relationship of Equation (1.17), we also store £ on the grid. 
Since s is independent of time, we can compute £ 1 prior to the loop for n . 
Additionally, in order to compute £ 1 , we solve Equation (1.4) for 0(z) and y/(z) for 
given boundary conditions using MATLAB’s bvp4c algorithm. 


By contrast, the magnetic field components are computed to coincide with points 
half way between grid points in both time and space, further referred to as the half-grid as 
seen in Figure 6. Again, this allows for the method to be second order but only requires 
storage of one field component per grid cell. 


H:Zl2=HM + U2,t„ +m ) 


H n y Zl2=H y (z k+ll2 ,t n+l/2 ) 


(1.25) 

(1.26) 


Both the electric field and magnetic field are initialized to zero throughout the entire 
computational domain. 
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Figure 6. Staggering of grid and half-grid. 


We now present the method for maintaining the update scheme current for the 
value of D x . We use a second order approximation of Equation (1.14) at t = t n+y2 and 

~ — Zk > 


8D 


dt 


dH 


= —c n 


z=z k 


dz 


z=z k 


D 


x.k 


-Kk 


A t 


+ 0 (aT ) — — Cq - 


H 


n+ 1/2 
y,4+1/2 


-H 


n+ 1/2 
y,4-1/2 


+ 0(az 2 ). 


AZ 


The values of D" k , H" )'^ /2 , and H" + k _ y2 are known from the previous time step, 
therefore we can solve for D at next time step in order to get an update scheme that is 
second order in both time and space, 


£>" t 1 = D" 

X.K X, 


At / u 

c o (H 

A Z 


n+ 1/2 
y,k+ 1/2 ' 


Zjn+ 1/2 

' 17 y,Jfe-l/2 


) + 0(At 2 ) + 0(AZ 2 ). 


(1.27) 


Similarly for Z), 


= 3U- 

/Tjn+1/2 Tj n+\/2 \ 

~ C 0 V 77 *,*+1/2 77 1/2 / ‘ 

(1.28) 


AZ 



= 

D n+i 

c xx,k^x,k 

+ +; 1 n n+1 

^ t’xyM y,k 

(1.29) 

= 

s { D n+X 

b xy,k 1J x,k 

+ £ ~ l D" +1 

^ t’yyk y,k 

(1.30) 
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From Equation (1.17), 













where z k ' is the inverse of z at z k . Note, however, that E" + k = ' k D" + k ' + a\', k D" + k is 
nonzero, although it not required to compute except for when it is necessary to compute 

|Er 1 - 


The update scheme for H is 


rr n+3/2 _ TTII+ 1/2 . ( Z7” +1 _ 

n x,k+V2 ~ n x,k+V2 ~ I ~ C 0 ' C ».W ' 

AZ 


rr n+3/2 Tjn+ 1/2 ^ t-,/?+ 1 Z7«+l\ 

AI 7,, r /t — Z3 yjA , +1/2 C 0 x^x.k+l ~ +.// ) ' 


1 y,*+l/2 


AZ 


(1.31) 

(1.32) 


Slight modifications to this update scheme must be made outside the liquid crystal 
layer. In the vacuum and PML, Equations (1.29) and (1.30) are simplified since z 1 is the 
identity so that E" k = D”. The update schemes for D" and H"^ within the PML are 
more complicated due to the addition of an absorbing medium. 


This scheme applies to both the TF region, where it updates the total field, and the 
SF region, where it updates the SF, but at the boundaries between the two regions, the 
scheme must be altered to avoid approximating the spatial derivatives with one total field 
value and one scattered field value. To update D" + k at the bottom transition, we need 

the total field value H n+ , 112 , but since this is in the SF region, it represents the value of 

y> K scatbot- 1/2 G 1 

the scattered field. Therefore we must adjust the equation using the known incident field 
as 


jyt i+l _ jyi ^ / jj n 1 - 1/2 jj n + 1/2 

x ^ scat hot x ’kscat hot 0 y ^scatbot+l/2 ^ ^scatbot-112 itlC,y ,k sca tbot-l/2 

A Z 


h:;' 1 :, ) 


D /i+l rv7 . ^ / rjn+ 1/2 r_r77+1/2 rr77+1/2 \ 

v k AZ k Cn (AZ r k AZ x k AZ ■ / ) . 

y~ K scatbot y~ K scatbot u . — x ’ K scatbot+l/2 x ’ K scatbot-V2 '' ' ’ x ’ K scatbot-V2 


AZ 


The magnetic field is similar with the exception that the electric displacement field value 
that is updated is in the TF region whereas the magnetic field value that is updated is in 
the SF region. This results in a magnetic field value update 


H 


n+3/2 


- 1/2 


_ rr 77+1/2 A f ( T?n+l 

*’Ka,,b„,- 1 / 2 ' rL o 

A Z 


-E 


!+l 

k „„ 


-1-^7 


n +1 
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jjn+312 _ rrzi+1/2 _ (E n+ ^ _ E n+ * _ E n+ ^ 'j 

ykscatbot-V 1 yKatbot-U 2 0 y^scatbot y^scatbot — ^ * nc scatbot 

A Z 

This is similarly done at the top transition. 

In order to truncate the computational domain without introducing nonphysical 
reflections, we must introduce the PML in order to absorb the wave as it leaves the 
computational domain. The absorption properties of the layer are characterized by the 
absorption coefficient cr. The behavior of a and its impact on the portion of the beam of 
interest are such that we absorb the beam slow enough so that the PML does not act as a 
hard boundary and reflect the beam. Conversely, it must absorb the beam quick enough to 
minimize the magnitude of reflection within the layer and in turn cause any unwanted 
reflections. 


Choosing the absorption characteristics of cr allows us to imitate endpoints of the 
computational domain that go on to infinity. 


a(x) = <j B 




V M J 


where x is the depth into the PML and d is the width of the PML. 

Again, we present the method for maintaining the update scheme current for the 
value of D x , this time however, inside the PML. We use a second order approximation of 

Equation (1.14) at t = t n+V2 and z-z k . 


dD x _ 
dt 


+ crD = —c n 


y 

dz 


At 


+ CT, 


DZ'+D" 


= ~c, 


jt 71+1/2 jjr 71+1/2 
** y,k+V2 ~ Jrl y,k- 1/2 

AZ 


j^n+1 _ 2 A t(Ji 


2 c n At 


2 +a7ct. 


k jyn sjjn+ 1/2 _ jjn+112 


k n, , . v,k+l/2 y,k-U2 

Z+At(T k A z 


). 


(1.33) 


Similarly for D , H x , and H 

2 -At a, 


K = 


2 +a7ct, 


2c 

k _ 0 

^ x,k 


A t 


2+A t(J, A z 


(H 


ti+1/2 rjn+l/2 \ 
x,k+ 1/2 ~ n x,k-\!2 > 


(1.34) 
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Tjn+m _ 2 Atcx k+\I2 tj»+1/2 
n x,k+ 1/2 — n x,k+ 1/2 


2C„ 


At 


2 +Afcr, 


*+ 1/2 


IT n+m _ 2 Atcr k+ll2 IT n+ 1/2 
v,1+1/2 0 77 y,1+1/2 ' 


2+ A ^CT lt+1/2 AZ 

2c n 


(^,l+ 1 -^,l) 


A? 


2 +AZcr, 


1 + 1/2 


2 +A/C jt+1/2 A z 


(Kti-Kt)- 


(1.35) 

(1.36) 


The incident beam at z start = k start /\z is introduced in the form of 

E" = E 1 sin(-(ynAt) 

X ’ K start 

E l 1 JB „ =0 


(1.37) 


Time, t, depends on the ^ of the material. This wait allows for initialization 
transient to pass and allow beam to travel one complete period such that intensity. 



CO 

i= J |exh|* can be calculated. 

t 


In order to calculate the intensity of the transmitted beam, we must first calculate 
the components of the Poynting vector S. Knowing that H_ - 0 



This requires that the components be known at the same space and same time. By the 
very nature of our staggered grid set up, this requires us to interpolate between grid points 
or half-grid points for E or H as necessary to compute the value of S. We choose to 
compute S on the grid so we must compute H to correlate with the space and time of the 
grid from the half grid. For example in the z- component, 


5 


n 

Z,k 



where H is the recomputed H from the half-grid to the grid and 


H " _ 1 / A ( tj+1+1/2 . tt w—1/2 iin+1/2 iin-1/2 \ 

k - 1/4 ( H i+l/2 +±1 *+l/2 +±1 *-l/2 +±1 *-1/ 2 J- 


Once S is known, we then use the trapezoidal method to approximate the intensity 
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where 7J represents the start time of the last period. This calculation is computed for a 
value of k that satisfies k glasstop+1 <k< k scattop+{ . 
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II. ANALYTICAL SOLUTION AND NUMERICAL VALIDATION 


When the same anchoring conditions are used on the top and bottom plates, we 
refer to Equations (1.4) because these equations have constant solutions throughout the 
liquid crystal layer. This is to say that when the anchoring conditions are matched on the 
top and bottom plates, this results in molecule orientations that are homogeneous 
throughout the LC layer. In this scenario, (// and 0 are constant across the LC layer. 

In this case, the electromagnetic problem can be solved analytically. In this 
chapter, we present a derivation of this solution and then we compare the results to the 
numerical solution in order to determine the optimal grid cell size for the given h . The 
basis for determining the optimal grid cell size is to ensure that we solve the problem 
with reasonable accuracy but also within a reasonable computational time frame. 

A. ANALYTICAL SOLUTION 

In order to solve the problem analytically, we introduce the known incident beam, 
with components E 7 and H 1 , from Equation (1.37), and four unknown resultant beams 
that we will solve for in terms of the incident beam. The reflected beam, E A ’ and H s , 
exists on the same region as the incident beam but travels in the opposite direction. The 
transmitted beam, E 7 and H r , exists in the opposite region as the incident beam and 
travels in the same direction. Internal to the LC layer are two beams. The first is 
comprised of E L+ and H i+ and travels in the positive z -direction whereas the other 
beam is comprised of E 7 ' and H 7 and travels in the negative z -direction within the 
layer. These are depicted in Figure 7. 
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Figure 7. Analytical beam components. 


In this section, we will use the complex forms of the electric and magnetic fields 
and we also introduce a phase shift that allows us to match these results with the results 
of Equation (1.37). Therefore, our fields take the form 

Ee _i£ffl g~ i( - k, > Zs ’ an+ ’ r/ ^ 


Plugging these forms into Equations (1.11) through (1.13) gives us the following 

-ifl*-E = c 0 VxH (2.1) 

ia>H = c 0 VxE. (2.2) 


From Equation (1.37), the incident beam is 


v ° 

1_ 

e' v , 

H 7 = 

o faf 

1_ 

1 

o 

1_ 



1 

o 

_1 


(2.3) 


where k 0 = ci)Je 0 ju 0 = — is the wave number in vacuum. The transmitted and reflected 

c o 

beams take the forms 
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'Ey' 

E x = 

K 

e~ ik ° z . 

H x = 

~E* 


0 



0 





'-El' 

E r = 

K 

e ik ° z . 

H r = 

El 


0 



0 


(2.4) 


(2.5) 


To express the unknown beam within the layer, we must first determine the wave 
numbers and the basis vectors for waves supported by the anisotropic LC medium. We 
need to define this constant vector E = E 0 e' fe and the corresponding H = H 0 e' fe for the 

unknown constant vectors, E 0 and H 0 and the wave vector k = (0,0, A:) 7 where k is an 

unknown wave number that is either positive or negative. Using these forms, Equation 
(2.1) and Equation (2.2), are now 


coz • E 0 = c 0 k x H 0 

(2.6) 

ruH 0 = c 0 k x E 0 

(2.7) 


We now solve Equation (2.7) in terms of H 0 and substitute it into Equation (2.6). 

Q 

-ox,- E 0 =c 0 kx(—kxE 0 ) 
co 

2 

£-E 0 =kx(kxE 0 ) =k(k-E 0 )-(k-k)E 0 
c o 

= kk-E () -k 2 E 0 

, co 2 

(k 2 I-kk-—-£)-E 0 = 0. (2.8) 

c o 

co 2 

The 3x3 matrix, k 2 l -kk- -e , must be singular, therefore we can set its determinant to 

c o 

zero and solve for k . 

2 

det(k 2 I - kk -——z) = 0 
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£ J* —rK £ B + f w) £ a -4 - £ i^ 2 


2 Ly xx yy / zz xy yz J 


H --[(£ 6 * -£ 2 )s -£ £ 2 -£ £ 2 +2s £ £ 1 = 0 

4 L v xx yy xy' zz xx yz yy xz xy yz xz J 


This is a quadratic equation in k 2 . Recalling that £ = (s i +A£nn) and knowing that nis 
a unit vector, we arrive at 

6 xx = £ ± +*£n; 

£ yy = £ L +A£n l 
£„ = £ ^ ~\~A£U_ 

£ xy = A £ n x n y 
£ xz = A£ n x n z 
£ yz = A£ n y n z 


Solving with substitution, we obtain 


2 (O 2 e l +£ 1 A£(l + n:)±£ 1 A£(l-n:) 


k - —— 

~ c 2 
L o 


2(e, +a £n 2 ) 


Utilizing 2s 2 + 2 e l as = 2s ± (s +as) and a£ = £,, -£ L , we get 


CO 


C0 2 £ l £, 


C0 2 £ i £„ 


cl(£ ± +A£n 2 ) cl{£, +A^sin 2 0 ) cl{£, cos 2 0 + £,, sin 2 d) 


This results in two distinct wave numbers that the wave can propagate through the 


LC layer. We stall with the ordinary wave number, k = — = — n , where n = ^[c~ 


is the ordinary index of refraction. The corresponding eigenvector of Equation (2.8) gives 


us a basis vector for the electric field 


E 0 = E ± = cos y/ 


The corresponding basis vector for the magnetic field is 
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H 0 =H ± = —kxE 0 =-n ± 
co 


cosy/ 

sin^ 

0 


The extraordinary wave number is k g - — - 1 11 --—=— n 0 , where 

c 0 V s ± cos' 6 + s n sin" 6 c Q 

£ | £u 

n a = ------— is the extraordinary index of refraction. In the case that 

y s L cos 2 0 + s^ sin 2 6 


6 = —, n e = n and so the material is effectively isotropic for this anchoring condition. 
We have the corresponding basis vector for the electric field 


E (I =E, 


cos y/ 
sin y/ 

-sin^cos^ 


H() H g Tig 


-sin^ 
cos y/ 
0 


These basis vectors allow us to express the electric and magnetic fields within the layer 
as 

e l± = E^ E±e ^z + E L± E ^ e ±ik e z 

H l± = ±E[ ± n L e ±ik ^ ± Ef YL d e ±ik <> z . 

We now have the general solution for the electromagnetic field but we must now 
apply the boundary conditions at the two interfaces at z — 0 and z = h. The boundary 
conditions come from forcing the continuity of the tangential components at the 
boundaries. That is that E x ,E y ,H x , and H must be continuous at z =0 and z-h . This 

creates a system of eight equations that will be used to determine the eight unknowns of 
those equations. 
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We will then look at the equations at z - 0. 

E{ + E r x = - sin y/{E[ + + E L ~ ) + cos y{E L ; + E L e ~ ) (2.9) 

E R = cosy/(E k+ + E^) + sm\//(E g + + E g ~) (2.10) 

E R = n L cosy/(-E k+ + E l f) + n 0 sin ^(-E g + + E'f) (2.11) 

E' x -E r x =n L siny/(-El + +E L L -) + n d cosy/(Ef-E l 9 ~). (2.12) 


And similarly, we will look at z = h 

E T x e ik ° h = -sin i//(E^e ikJl + E L ~e ~ ikih ) + cos i//(E L g + e ikeh + E^e'* 6 

E T y e ik ° h = cosi//(E k+ e ,k h + E l fe" k ' h ) + sim//(E I g + e ikeh + E g ~e^ h 

-E T y e ik ° h = n cos i//(-E k+ e iklh + E'fe~ ik ' h ) 

+ n 0 sin y/(-E L e + e ikeh + E^e~ ik ° h ) 

E T x e ik ° h = n ± sin if/(-E^e ,k ' h + E^e lkh ) 

+ n g cos i//(E k+ e lVl - £''V * 9 *) 


") 

) 


(2.13) 

(2.14) 

(2.15) 


(2.16) 


We will meticulously walk through the derivation of an equation through the 
summation and subtraction of Equations (2.9) through (2.16). This will allow us to 
simplify the system of equations to four equations with four unknowns. 


Equation (2.13) plus Equation (2.16) results in 

2 E T X e' k " h = sin i//\-E ,A e ,k h (1 +nj- E k e ,k h (1 - n ± )] 
+ cos y/[E g + e ikeh (1 +n g ) + E^e~ ikt,h (1 - n g )] 


(2.17) 


Equation (2.13) minus Equation (2.16) results in 

0 = - sin y/[E L 'e' k ' h (1-n ) + E L e ,k ' h (1 + n L )] 

+ cos y/[E L g + e iksh (1 -n g ) + E^e'** (1 + n g )] 

Equation (2.14) minus Equation (2.15) results in 

2 E] e ik ° h = cos i//[E' ± + e ik ' h (1 + n ± ) + (1 - n ± )] 

+ sin i//[E g + e' kgh (1 + n e ) + E L e ~e~ ikeh (1 - n g )] 


Equation (2.14) plus Equation (2.15) results in 
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(2.19) 



0 = cos y/[E''e‘ k h (1 -nj + E l [e lk ' h (1 + n )] 
+ sin ys[E k+ e ikeh (1 -n e ) + E'fe^ (1 + n e )] 


( 2 . 20 ) 


Equation (2.20) times cos y/ plus Equation (2.18) times - sin yr results in 

0 = sin 2 y/\E['e' k ' h (1 -nj + E L e ,k h (1 + n L )] - cos y/ sin y/[E L g + e ik ° h (1 -n g ) + E k e ik ° h (1 + n g )\ 
+ cos 2 y/\E k 'e' kh (1-/7) + E l e ,k h (1 + n ± )] + cos y /sin y/[E k+ e ik(,h (1 -n g ) + E k e ,k " h (1 + n g )] 
= E k+ e ,k ' h (1 -nj + E'fe~ ,k ' h (1 + nj 




e~ /2 ^(l + nj 
(1 -nj 


( 2 . 21 ) 


Similarly, 


E k+ = -E a 


- i2k ° h (l + n g ) 

(1 -n g ) 


( 2 . 22 ) 


Equation (2.9) plus Equation (2.12) results in 

2 E[ = -sini//[E k+ (1 + n ± ) + E k (l- n ± )] 

+ cos y/[E g + (l + n g ) + E g (l-n g )\ 

Equation (2.10) minus Equation (2.11) results in 

0 = cos if/[E k+ e lk ' h (1 + nj + E'fe~ ,klh (1 - n L )] 
+ sin y/[E L g + e ik ° h (1 + n g ) + E k ~e^ h (1 - n g )] 


(2.23) 


(2.24) 


Equation (2.23) times cos y/ plus Equation (2.24) times sin ^/results in 

2 E' x cos y/ = E L e + (1 + n g ) + E k ~(1 -n g ). (2.25) 


Equation (2.24) times cos y/ minus Equation (2.23) times sin ^/results in 

-2 E[ sin y/ = E k+ (1 + n ± ) + E[~ (1 - n ). (2.26) 

Substituting Equation (2.21) into Equation (2.26) results in 

-2 E[ ssn ¥ = El-m-n,)- - ( "^ ]. 

(1 ~nj 

This leads us to a solution for the electric field related to the ordinary index of 
refraction within the liquid crystal layer in both the positive and negative z -directions 
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2 E[ sin^(l-n ± ) 
-nj -(1 + nj e 1 

2 E{ sin i//(l + n ± ) 

a \2 ilk \h /i i \ 2 

-nj e -(l + n ± ) 


Similaiiy for the extraordinary component 


E l ~ = 

^9 


2 E[ cos^(1 -n g ) 
(l-n e ) 2 -(\ + n g ) 2 e- i2k o h 


f l + _ ~2E x cosy/(l + n g ) 

6 (l-n d ) 2 e i2k ° h -(l + n 9 ) 2 


We can now get the transmitted electric field in both the x- and y- directions by 


substituting E g , E' g + , E[ , E 2+ into Equation (2.13) 


-± ’ 
pL+ ik | h 


- sin y/(E' ± + e ,K,n + E k ~e~ ik±h ) + cos y/(E' e + e ,Ken + E'fe '**) 


L+ iknh 


■pL— —ikah\ 


E = 


iknh 


sin y/ 


e ik ° h 


2 E’ x sin ^[(1 + n x ) - (1 - n ± )] 
(1-/7, )V M -(1 + /7, )V' M 


+ 


cos^ 


Aknh 


2 ^.v cos y^[(l — n g ) — (1 + 77 ^)] 
(\-n g ) 2 e ik ° h -(\ + n e ) 2 e- ik ° h 

n L E[ sin 2 y/ 


iknh 


(1 -n.) 2 e ,k±h ~(\ + n,ye 


2 -ik.h 


n g E{ cos 2 y/ 


(1 - n g ) 2 e lk(,h -(\ + n g Ye 


2 -ik n h 


El = 


ae: 


x iknh 


n ± sin y/ 


n e cos y/ 


/-i . \2 -ik\h \2 ik\h /i , \2 -iknh /-t \2 u 

(l + n ± ) e 1 -(1 -n ± ) e 1 (1 + n g ) e e -(1 -n g ) e 


Similarly, 


El = 


AE. 


y iknh 


n g cos y/ sin y/ 


n ± cos y/ sin y/ 


(1 + 77,)V' V ‘ -(1 -n,)V v (1 + nXe" K * -(1 -nXe 


\2 ikgh 


2 —ik i h 


\2 ik.h 


We now have solutions for the electric field of the transmitted beam with respect 
to the incident beam which is a known quantity. Our next step is to derive the intensity of 
the transmitted beam with respect to known quantities. 
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1 CO 

I - ——-— J \ExH\dt 


i T = - XE T .E T *=- ^ 4 e t e t *+e t e t *) 

2^ 0 - - 2 v// o ' ' y v 


r 16 / £* j 2 

' "it ‘ 


. 2 

n L sin y/ 


Tin cos \j/ 


/i \ 2 —ik « /•« -.2 ikh -.2 ~tk.h \2 ikJi 

[_(1 + Hj_) e -(l-rc ± ) e (l + rc^) £ -(1-n^) e 

. 2 2 
n_L sin y/ n,Q cos y/ 

(1 + n, ) e — (1 —w, ) e (H-tu) e -(1 — £ 


sin y/ cos y/ 


n L cos y/ sin y/ 


/I i \2 —ik h \ 2 ifc./i /I ■ \2 ik h ✓ -« \ 2 ifcJ 

(1 + ^) e — (1 — ) e (1 + w , ) e — (1 - w , ) e * 


n e sin y/ cos y/ 


n L cos y/ sin y/ 


= 8\—(e[) 2 


/*1 \ 2 JK.rt \ 2 ik h \ 2 (fc h \ 2 /K 

(1 + tin ) e - (1 -tin) e (1 + n , ) e -(1 -n,)e 


n L sin y/ 

(1 + n ± ) 4 - (1 + n ± ) 2 (1 - n ± ) 2 (e i2 ‘ kJ ' + e^") + (1 - n ± ) 4 


cos i// 


0- + n g ) 2 -Q- + n„) 2 (l-n e ) 2 (e 6 + e 6 ) + (1 - n 9 ) 4 


2.2 22 
s , , n sin w n„ cos y/ 

_4 — —(E 1 y --- - -+- e - - - - 

V H o ' \ + 6n 2 +n ± 4 — (l-« ± ) 2 cos(2& ± /i) l + 6« e 2 +n g 4 — (l-« g 2 ) 2 cos(2k e h) 


2.2 22 

, «_l sin y/ n e cos y/ 

— ol --- - --- 1 ---—-- 

1 + 6 n ± ~ +n L — (1 -n ± ) 2 COS(2 k h) l + 6n ff " +n 4 — (!-«')' COs(2 k h) 


where the intensity of the incident beam is I =— —(E 2 ) 2 . Solving for the remaining 

- V /'„ 

unknowns with respect to known values and summarizing those already solved for, we 


/' = 8 / 


n ± sin“ y/ 

4 „ \ 


1 + 6 n L +n L -(l-n ± ) cos(2k L h) 
n g 2 cos 2 (// 

1 + 6n g 2 + n g 4 - (1 - n 2 ) 2 cos(2 k g h) 


(2.27) 
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(2.28) 


El = 


4 El 


iknh 

e 0 


n sin i// 


(\+nJe- ,k ' h -(\-n L fe u 


n g cos y/ 


(1 + n g ) 2 e~ ikeh -(1 - n 0 f e ,k(,h 


El = 


4 El 


y iknh 


n 0 cos y/ sin y/ 


(l + n d ) 2 e- ik ° h -(l-n e Ye 
n L cos y/ sin y/ 

/i , \2 -ik,h \2 ik\h 

(1 + 77,) e L -(1-77,) e 1 


2 ik„h 


Efr=- 


2E\ sin^l-Tj^) 

(1-/J ± ) 2 -(1 + /1_ L )V'' 2M 


e l : = 


2E\ sin^(l + 77 ± ) 

(1-/J ± )V 2M -(1 + /1 ± ) 2 


?l+ _ -2£ A / cosy7(l + /7 g ) 

(l-77,) 2 ^-(l + /7,) 2 

2E\ COS^(l-77 0 ) 


e l ~ = 


(l-77,r-(l + 77,)^ 


2 —ilkah 


E R = 


(1 -n e 2 )E' x cos 2 y/(l -e~ 2lkr>h ) 


(l-n s y-(l + n g y e 


2 -2 iknh 


(1 - n 1 z )E‘ t sin" y/{\ - e i2k±h ) 

(l-w 1 )V 2 *- L ‘-(l + /i ± ) 2 


cos y/ sin i// 


(\-n e 2 )(\-e~ i2k ° h ) 

(\-n 9 y-(\ + n g ) 2 e- 2ik ° h 

, (l-/i ± 2 )(l-g f2M ) 

(l-77 1 ) 2 e' 2M -(l + 77 1 ) 2 


(2.29) 

(2.30) 

(2.31) 

(2.32) 

(2.33) 

(2.34) 

(2.35) 


Equation (2.27) gives an explicit formula for the intensity of the transmitted beam. We 
will use it to validate our FDTD codes in the next section. 
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B. 


VALIDATION OF NUMERICAL CODE 


Before we can compare the analytical and numerical solutions, we must first 
determine, through the FDTD method, the optimal mesh size that takes into account the 
relationship between our time step and our spatial step that will give accurate results. We 
compute multiple transmitted intensities with a varying mesh grid size and show them in 
Figure 8. This illustrates that all solutions converge to a mesh grid with a value of grid 
points approximately N=3 100. It is through this that we can say with confidence that the 
3200 grid points will provide a solution that is consistently accurate. However, we must 
also look at acceptable error and the time required computing the numerical solution. 


Intensity vs Grid Size with Psi= 0 degrees 



Figure 8. Transmitted intensity as compared to grid spacing with y/ — 0. 

Additionally, we look at the overall error between the known true solution and the 
numerical solution for different grid mesh sizes. Table 1 is a collection of various values 
of intensity errors between the true value and the numerical value for the given N, y/, 
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and where 6 varies from 0 to ir 12 in increments of n / 20. We have provided the 
maximum, minimum and average values for each of the values of N. Additionally, we 
have provided the percentage error which is close to the value of the difference in 
intensity. 
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Mesh to True Solution Error 

Mesh to True Solution Percent Error 



200 

400 

800 

1600 

3200 

200 

400 

800 

1600 

3200 

(0,0) 

MIN 

0.0232 

0.0037 

0.0052 

0.0012 

0.0004 

2.49% 

0.47% 

0.52% 

0.12% 

0.04% 

MAX 

0.2667 

0.0797 

0.0294 

0.0135 

0.0064 

26.67% 

8.01% 

3.19% 

1.40% 

0.67% 


AVG 

0.1408 

0.0479 

0.0183 

0.0077 

0.0035 

15.18% 

5.15% 

1.98% 

0.83% 

0.38% 

(n/20 ,0) 

MIN 

0.0259 

0.0043 

0.0052 

0.0012 

0.0004 

2.77% 

0.53% 

0.52% 

0.12% 

0.04% 


MAX 

0.2635 

0.0784 

0.0288 

0.0131 

0.0063 

26.35% 

7.88% 

3.12% 

1.36% 

0.65% 

AVG 

0.1389 

0.0470 

0.0179 

0.0075 

0.0034 

14.92% 

5.03% 

1.93% 

0.81% 

0.37% 

(/r/10,0) 

MIN 

0.0338 

0.0060 

0.0052 

0.0012 

0.0004 

3.38% 

0.60% 

0.52% 

0.12% 

0.04% 

MAX 

0.2543 

0.0747 

0.0271 

0.0121 

0.0058 

25.43% 

7.47% 

2.71% 

1.21% 

0.58% 


AVG 

0.1334 

0.0441 

0.0167 

0.0070 

0.0032 

13.34% 

4.41% 

1.67% 

0.70% 

0.32% 

{in /2O,0) 

MIN 

0.0375 

0.0086 

0.0052 

0.0012 

0.0004 

3.38% 

0.60% 

0.52% 

0.12% 

0.04% 


MAX 

0.2399 

0.0690 

0.0245 

0.0105 

0.0050 

25.43% 

7.47% 

2.71% 

1.21% 

0.58% 

AVG 

0.1247 

0.0397 

0.0149 

0.0062 

0.0028 

13.34% 

4.41% 

1.67% 

0.70% 

0.32% 

Or/5,0) 

MIN 

0.0075 

0.0096 

0.0052 

0.0012 

0.0004 

0.75% 

0.96% 

0.52% 

0.12% 

0.04% 

MAX 

0.2218 

0.0618 

0.0211 

0.0089 

0.0041 

22.18% 

6.18% 

2.11% 

0.89% 

0.41% 


AVG 

0.1137 

0.0341 

0.0125 

0.0051 

0.0023 

11.37% 

3.41% 

1.25% 

0.51% 

0.23% 

(*/4,0) 

MIN 

0.0128 

0.0010 

0.0040 

0.0012 

0.0004 

1.28% 

0.10% 

0.40% 

0.12% 

0.04% 


MAX 

0.2016 

0.0538 

0.0174 

0.0071 

0.0032 

20.16% 

5.38% 

1.74% 

0.71% 

0.32% 

AVG 

0.1085 

0.0279 

0.0100 

0.0040 

0.0018 

10.85% 

2.79% 

1.00% 

0.40% 

0.18% 

(3^/10,0) 

MIN 

0.0190 

0.0002 

0.0012 

0.0008 

0.0004 

1.90% 

0.02% 

0.12% 

0.08% 

0.04% 

MAX 

0.1814 

0.0457 

0.0137 

0.0053 

0.0024 

18.14% 

4.57% 

1.37% 

0.53% 

0.24% 


AVG 

0.1126 

0.0231 

0.0074 

0.0029 

0.0013 

11.26% 

2.31% 

0.74% 

0.29% 

0.13% 

{In! 20,0) 

MIN 

0.0651 

0.0102 

0.0007 

0.0000 

0.0001 

6.51% 

1.02% 

0.07% 

0.00% 

0.01% 


MAX 

0.1631 

0.0384 

0.0103 

0.0036 

0.0016 

16.31% 

3.84% 

1.03% 

0.36% 

0.16% 

AVG 

0.1216 

0.0246 

0.0052 

0.0019 

0.0008 

12.16% 

2.46% 

0.52% 

0.19% 

0.08% 

(2^/5,0) 

MIN 

0.1020 

0.0192 

0.0019 

0.0000 

0.0001 

10.20% 

1.92% 

0.19% 

0.00% 

0.01% 

MAX 

0.1484 

0.0326 

0.0076 

0.0023 

0.0010 

14.84% 

3.26% 

0.76% 

0.23% 

0.10% 


AVG 

0.1288 

0.0261 

0.0049 

0.0012 

0.0005 

12.88% 

2.61% 

0.49% 

0.12% 

0.05% 

(9x/2O,0) 

MIN 

0.1262 

0.0251 

0.0043 

0.0008 

0.0002 

12.62% 

2.51% 

0.43% 

0.08% 

0.02% 


MAX 

0.1388 

0.0288 

0.0058 

0.0015 

0.0005 

13.88% 

2.88% 

0.58% 

0.15% 

0.05% 

AVG 

0.1334 

0.0270 

0.0051 

0.0011 

0.0004 

13.34% 

2.70% 

0.51% 

0.11% 

0.04% 

(*/2,0) 

MIN 

0.1352 

0.0274 

0.0052 

0.0011 

0.0004 

13.52% 

2.74% 

0.52% 

0.11% 

0.04% 

MAX 

0.1352 

0.0274 

0.0052 

0.0012 

0.0004 

13.52% 

2.74% 

0.52% 

0.12% 

0.04% 


AVG 

0.1352 

0.0274 

0.0052 

0.0012 

0.0004 

13.52% 

2.74% 

0.52% 

0.12% 

0.04% 


Table 1. Absolute and percent error between the true intensity and the numerical 

intensity for the given N, (//, and 0 . 
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The information in Table 1 can be summarized in Table 2 such that a more 
reasonable comparison can be done for which value of N to choose. In Table 2, we have 
averaged all of the minimum, maximum and average values for the given N as well as list 
the maximum and minimum values for that value of N. 


Mesh to True Solution Error (ty/, 6 ) Mesh to True Solution Percent Error 



200 

400 

800 

1600 

3200 

200 

400 

800 

1600 

3200 

AVG MIN 

0.0511 

0.0101 

0.0035 

0.0008 

0.0003 

5 . 09 % 

0 . 99 % 

0 . 35 % 

0 . 08 % 

0 . 03 % 

AVG MAX 

0.1774 

0.0465 

0.0147 

0.0060 

0.0028 

17 . 87 % 

4 . 71 % 

1 . 52 % 

0 . 62 % 

0 . 29 % 

AVG AVG 

0.1139 

0.0293 

0.0091 

0.0035 

0.0015 

11 . 57 % 

3 . 00 % 

0 . 94 % 

0 . 36 % 

0 . 16 % 

MIN 

0.0075 

0.0002 

0.0007 

0.0000 

0.0001 

0 . 75 % 

0 . 02 % 

0 . 07 % 

0 . 00 % 

0 . 01 % 

MAX 

0.2667 

0.0797 

0.0294 

0.0135 

0.0064 

26 . 67 % 

8 . 01 % 

3 . 19 % 

1 . 40 % 

0 . 67 % 


Table 2. Summary of intensity error and intensity percent error from Table 1. 


The error introduced by selecting a value of N=800 is an acceptable error for purposes of 
evaluating the overall impact of varying our anchoring conditions. All figures, unless 
otherwise stated, are generated using a mesh grid with A=800. 

We have shown that we are able to compute the known true solution for the 
electromagnetic wave with matched anchoring conditions. We show through the FDTD 
method and computer computation time requirements that the most efficient mesh grid 
spacing for accuracy to the known true solution is iV=3200 but this is not the most 
efficient when compared to acceptable error and time computational requirements. 
Figure 9. depicts the comparison of the normalized intensity of both the numerical and 
analytical solution using a mesh grid with A=3200. This is achieved by setting i// -0 
and varying 9 from 0 to n / 2. It should also be noted that the graph depicts multiple 
values of 9 that produce a maximum intensity but only one value that produces a global 
minimum intensity. We will look at this closer later. 
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Normalized Intensity 


Matched Anchoring Transmitted Intensity with vj; = 0 



Figure 9. True versus numerical solution. 
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III. NUMERICAL SIMULATIONS 


A. MATCHED PSI 

As a point of convention, we will use the term intensity and transmitted intensity 
to mean the same intensity of the plane wave. Otherwise we will clarify the intensity of 
point. We start first by stepping away from the known solution of matched, y/ top = y/ bot 

and 0 = 9 hot , anchoring conditions and shift to a combination of mixed, y/ ^ y/ bot 
and/or 0 top 9 bot , and matched anchoring conditions on both the top and bottom of the 

LC layer. At this point we can no longer utilize the true analytical solution in our process 
as it is not valid for mixed anchoring conditions. We have set y/ — 0 and maintain that 
value for the duration of each individual data run. 

Figures 10 through 19 are similar in that the value of yr is held constant and each 
individual curve on the specific graph is represented by a value of 9 top ranging from 0 to 

n 1 2 in increments of n / 20. While holding 9 top at that value, we vary 9 bot from 0 to 

n 12 in increments of n /100. The important point to take away from the curve is not 
the data itself but more as to what happens to the shape and relative relationship of the 
curves as y/ is changed while maintaining matched conditions. 

It is from Figure 10 through Figure 19 that we see that the relative relationship 
between the curves on one graph do not change as compared to those relative 
relationships on another graph. What we see is the effective compression of the curves 
and the minimum intensity going up as y/ increases. We also note that the maximum 
value of transmitted intensity does not change. It is seen that as the value of y/ 
approaches y/ — n 1 2, the intensity solution converges to a value of one. This means that 
the layer is essentially transparent to the beam with that specific polarization. 
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Figure 10. Transmitted intensity - y/ matched with 6 bot ranging from 

0 to n! 2 in n! 20 increments. 


4i=rc/20 



Figure 11. Transmitted intensity - y/ matched with 6 bot ranging from 

0 to n! 2 in n 1 20 increments. 
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M<=7r/10 



Figure 12. Transmitted intensity - y/ matched with 6 hot ranging from 

0 to 7i 1 2 in 7i 1 20 increments. 


M <= 3 r /20 



Figure 13. Transmitted intensity - y/ matched with 6 hot ranging from 

0 to n I 2 in n I 20 increments. 
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\y=n/5 



Figure 14. Transmitted intensity - y/ matched with 6 bot ranging from 

0 to n 1 2 in n I 20 increments. 


4 >= t :/4 



Figure 15. Transmitted intensity - y/ matched with 0 bot ranging from 

0 to k! 2 in n 1 20 increments. 
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vji = 3 t ^10 



Figure 16. Transmitted intensity - t// matched with O hot ranging from 

to n! 2 in n / 20 increments. 


> t «= 7*/20 



Figure 17. Transmitted intensity - y/ matched with 6 hot ranging from 

0 to n! 2 in n! 20 increments. 
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»| i = 2 n /5 



Figure 18. Transmitted intensity - y/ matched with 6 bot ranging from 

0 to n 1 2 in n 1 20 increments. 


4J=9r./20 



Figure 19. Transmitted intensity - y/ matched with 6 hot ranging from 

to n / 2 in n! 20 increments. 
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Figure 20. Transmitted intensity - t// matched with 6 bot ranging from 

0 to k! 2 in k! 20 increments. 

This is a significant finding in that anytime we have y/ - n / 2, as shown in 
Figure 20, the intensity solution is independent of 6 . 

We now move forward and investigate the effects of setting the value of 9 to a 
fixed value, vary 0 bot from 0 to n / 2 in fifty increments and vary the matched value of 
y/. The results as shown in Figure 21 through Figure 31. show this relationship with 
each graph representing a different value of 9 top . The curve of y/ - 7t 1 2 is consistently 

at the top represented by the dashed blue line. Additionally the curve of y/ — 0 is 
consistently at the bottom of the collection of curves. These two curves represent the 
upper and lower bounds of the solution and this is consistent with the solution for 
intensity as shown in Equation (2.27). 
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Normalized Intensity Normalized Intensity 


Matched 0 to it/2 with otop=0 



Normalized transmitted intensity for y/ — 0 to y/ — n / 2 and 9 top fixed. 


Matched v 0 to rt/2 with otop=it/20 



Normalized transmitted intensity for y/ — 0 to yr — n / 2 and 0 top fixed. 
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Figure 22. 








Normalized Intensity Normalized Intensity 


Matched vy 0 to it/2 with otop=it/10 



Normalized transmitted intensity for y/ — 0 to y/ — n / 2 and 0 top fixed. 


Matched v 0 to */2 with otop=3»t/20 



Normalized transmitted intensity for y/ = 0 to y/ — n / 2 and fixed. 
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Figure 24. 











Normalized Intensity Normalized Intensity 


Matched ^ 0 to s/2 with etop=s/5 



Normalized transmitted intensity for y/ = 0 to iff = n / 2 and 0 top fixed. 


Matched ip 0 to s/2 with otop=s/4 



Normalized transmitted intensity for y/ — 0 to yr — n / 2 and 0 top fixed. 
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Figure 26. 










Matched 0 to e/ 2 with otop=3ir/10 



Figure 27. Normalized transmitted intensity for y/ = 0 to y/ - n l 2 and Q top fixed. 


Matched v 0 to it/2 with otop=7jt/20 



Figure 28. Normalized transmitted intensity for y/ -0 to y/ — n I 2 and 9 top fixed. 
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Matched 0 to i/2 with otop=2s/5 



Figure 29. Normalized transmitted intensity for y/ = 0 to y/ - n I 2 and 0 top fixed. 


Matched 0 to jt/2 with utop=2r/5 



Figure 30. Normalized transmitted intensity for y/ - 0 to y/ — n / 2 and 0 Hiji fixed. 
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Matched n> 0 to rJ2 with otop=9-/20 



Figure 31. Normalized transmitted intensity for y/ - 0 to y/ ~ n I 2 and 0 top fixed. 

In a focus to reduce the intensity of the transmitted beam, we can therefore utilize 
the lower bound of the scenario and for matched y /, we see that we get the maximum 
intensity reduction at y/ — 0 for any value or combination of 6 as compared to any other 
matched value for y/ . This allows us to, is the case of matched ;// , to set its value to 0 
and remove it as a variable impacting intensity. This then means that the overall intensity 
reduction is dependent on the value of 6 . 

By plotting both a surface plot and contour plot as seen in Figure 32 and Figure 
33, we see that there is a global minimum at lower values of 6 . The graphs are utilizing 
the same data which consists of the calculated intensity with yr — 0 and 50 data points for 
both 6 top and 6 hot . There is also a curvature of the valleys at the individual minimums as 
well as curvature at the peaks of Figure 32. This curvature can also be seen in Figure 33. 
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Figure 32. Surface plot of transmitted intensity with y/ — 0 for mixed 0 . 


Line 1 on Figure 33, the diagonal line that connects the lower left corner and the 
upper right corner, represents matched anchoring for all parameters where y/ top = i// hol and 

0 lop = 0 hol . The diagonal, Line 2, is perpendicular to Line 1 and contains the points where 

the values of 9 top and Q hot are such that they average to the same value for all points on 

that line. This is true for all lines parallel to Line 2. The area of interest that corresponds 
to the global minimum seen in Figure 30 is shown by Line 3 on Figure 33 and correlates 
to the area where the average is approximately 0.35 for 6 lop and 6 bot . In Figure 34, we 

take a closer look at the area of interest by taking a slice of Figure 32, where 0 top = 0 bot . 
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Figure 33. Contour plot of transmitted intensity with y/ — 0 for mixed 6 . 
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1 


Matched Anchoring Conditions with = 0 



Figure 34. Transmitted intensity for matched 6 and matched y/ — 0. 

As expected, we see that the global minimum is 0.3456 radians. From both the 
surface and contour plots, we want to find all values of 6 that are on that line that 
averages to 0.3456 radians. Figure 35 shows the transmitted intensity for different values 
of 0 where the average of Q and 6 bot is 0.3456. 
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Figure 35. Impact on transmitted intensity for mixed 9 around an 

average 9- 0.3456. 

It is clear that the best reduction that we can achieve for a plane wave of known 
polarization is to use matched y/ and matched 9 where the average of 9 - 0.3456. 
Additionally, when utilizing mixed theta for a matched psi, the minimum transmitted 
intensity is dominated by selecting 9 lop and 9 bot such that the average is still 0.3456 and 
the difference is near zero. 

B. MATCHED THETA 

We now take a look at a matched 9 with a mixed i//. We first look at a twisted 
orientation of the LC molecules by setting yr and y/ bot perpendicular to one another and 

then to swap those anchoring conditions. The results can be seen in Figure 36 through 
Figure 44 for matched values of 9 - 0, 0.3456, k1 4. Figure 36 shows that there is little 
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difference between the curves with regard to reduction in intensity. This is consistent 
with the associated contour and surface plot. The greatest reduction is a twisted structure 
established by setting y/ -nil and y/ bot = 0. 


Intensity with 0=0 Intensity with o=0 



Figure 36. Mixed y/ with matched 6-0 
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Normalized Intensily 
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Figure 37. Contour plot of 6 - 0 with mixed i// . 


Intensity with 0=0 
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Figure 38. Surface plot of 6 - 0 with mixed ;// . 
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The results of 6 -n 1 4 are similar in that the maximum reduction occurs with 
y/ = n / 2 and y/ hot = 0 however is should be noted that the reduction is about ten 
percent greater than with 0 - 0. 


Intensity with e=x/4 Intensity with e=?:/4 



Figure 39. Mixed t// with matched 6— n i 4 
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Normalized Intensity 


Intensity Contour at 0=s/4 



n>bot 


Figure 40. Contour plot of 0 - n / 4 with mixed y/ . 


Intensity with 0=si4 
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Figure 41. Surface plot of 6 - n / 4 with mixed ;// . 
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Finally, Figures 42 through 44 give a greater band to the value of 9 top as long as 

9 hot is bounded to a value <0.5. It should be noted that the maximum reduction occurs at 

0, op < 1 with 9 hot < 0.1. This is consistent with our findings for 9 = 0.3456 and i// = 0 for 

the matched conditions discussed before. The results shown here also suggest that there is 
more complicated relationship between the transmitted intensity and psi value. 


Intensity with e=.3456 


Intensity with e=.3456 



w top 


bot 


Figure 42. Mixed t// with matched 9 = 0.3456 
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Normalized Intensity 


Intensity Contour at 0=0.3456 
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Figure 43. Contour plot of 6 - 0.3456 with mixed \j/ . 


Intensity with 9=0 3456 
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Figure 44. Surface plot of 6 - 0.3456 with mixed t// . 
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IV. CONCLUSIONS 


A. SUMMARY OF RESULTS 

We have looked at the propagation of a plane wave through a LC layer, paying 
special attention to how the orientation of molecules within the layer affects transmission 
of an incident beam. Knowing that the anchoring conditions are what determines overall 
molecular orientation within the LC layer, our goal was to find the optimal anchoring 
conditions that minimize the transmitted intensity of a known incident beam through the 
LC layer. These four anchoring conditions, two anchoring conditions at each plate, are 
independent of one another. This eventually leads to an optimization problem in four¬ 
dimensional space. 

Given the complexity of this problem, we have constrained the scope of this study 
by limiting the incident beam to a specific beam of known parameters. We used an 
incident beam of known wavelength, polarization, perpendicular incident angle, and 
intensity. There exist other beams that may have a smaller or greater impact from the 
designed layer. The process of minimizing the intensity of a known beam is different 
than the process for the vast spectrum of possible beams in wavelength, polarization, and 
incident angle. This allows us to show preferential focus towards best reducing one 
specific beam intensity. 

We have used a FDTD code to propagate the plane wave through the 
computational domain. For orientations using matched anchoring conditions, we were 
able to validate that code by comparing the computed numerical solution to that of the 
known analytical solution. We validated the solution by comparing meshes consisting of 
200, 400, 800, 1600, and 3200 grid points, respectively. We looked at the error between 
the magnetic field solutions, the electric field solutions, and the intensity solutions 
through the entire computational domain. Knowing that the mesh grid with A=3200 
provided a more accurate solution, we chose to utilize a mesh grid with A=800 to 
conserve processing time without significant reduction in accuracy. 
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We first started looking at matched anchoring conditions for both y/ and 6 , of 
which the true solution is known through analytical means. We then looked at selecting 
fixed values of y/ and varying 6 . For yr — n / 2, the known true solution, for this beam, 
reduces to a solution that is not dependent on 6 . Consistent with the true solution, we 
observed that as matched yr approached x12, all solutions converged to maximum 
intensity. Ultimately for any matched value of yr , the solution is bounded by those 
solutions with yr - 0 and yr — n / 2. 

With both matched y/ and matched 0 , the relationship of intensity was consistent 
with that of Equation (2.27) and any value of matched y/ could be determined through 
the linear interpolation of the values of y/ — 0 and yr — n / 2 . Also, for all values of 
matched y/ , the value of yr = 0 was always the greatest intensity reduction regardless of 
the value of matched 0. This allowed us to set matched y/ to zero for both maximum 
intensity reduction and to remove it as a variable. We were then able to look at 0 
through the span of mixed and matched values and determine that the best value was an 
average 0 of 0.3456 with the a9 < .2 in order to show a maximum reduction of 
approximately 25 percent. This solution was verified by a run at a high mesh of 3200. 

This segued into the scenario of matched yr and mixed 0 , whereas mentioned 
before, the average of the thetas dominated the output. For matched 0 and mixed yr , our 
preliminary results show that there is more complicated relationship between the 
transmitted intensity and psi value. 

B. FUTURE/ONGOING WORK 

There continues to be studies in the specifics of this paper. Specifically to this 
study would be to expand to a 3-dimensional analysis of the problem. This would require 
a more efficient way to store data than was used in this study as memory was constantly 
an issue with certain data runs. Also, working toward a multilayer look would be a 
simple step forward by just treating this study as an individual layer. Multilayer also has 
the aspect of changing the thickness of the LC medium which was not done in our study. 
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Investigating into incident beams of different wavelengths and incident angles is 
very crucial to determine the feasibility of using a LC layer, or similar material, as a 
protective layer for DoD assets. Other areas of interest are coupled electromagnetic 
fields where the orientation can be adjusted as needed throughout the medium rather than 
just anchoring conditions at the boundaries. 
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APPENDIX 


A. MATLAB CODE FOR NUMERICAL SOLUTION 

This MATLAB code was used as the foundation in order to calculate and generate 
curves for the electric field, magnetic field, and intensity. The original code was created 
by Dr. Eric Choate and modified in order to suit our application. Many alterations were 
made for the various sections with this code as the baseline for computations. 

9 - 9 ' 2 - 9 ' 2 - 2 ' 2 - 9 ' 2 - 9 ' 2 - 9 ' 2 - 9 - 2 ' 9 ' 2 - 9 - 2 - 9 ' 2 - 9 ' 2 - 9 ' 2 - 9 ' 2 - 9 ' 2 - 2 ' 2 ' 9 ' 2 - 9 ' 2 - 9 ' 5 ' 9 ' 2 - 9 ' 2 - 9 ' 2 - 9 ' 2 - 2 ' 2 ' 9 ' 2 ' 9 - 2 - 9 ' 2 - 9 ' 2 ' 9 ' 2 - 9 ' 2 ' 9 - 2 - 9 - 2 - 9 - 2 - 9 - 2 - 9 - 9 - 9 - 2 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%This implements the Finite Difference Time Domain (FDTD) algorithm for 
propagation of the electric and magnetic fields of a plane wave across 
a system of liquid crystal polymer (LCP) domain. 

%The purpose of this code was to test the convergence of the algorithm. 
That is, by using a finer grid, we should get a better approximation of 
the true solution. FDTD is second order, and so we expect that the 
error should scale like dx A 2, where dx is the grid spacing. This code 
solves the same problem for six successively finer grids. On each 
grid, two different waves are propagated. The incident beam Exinc and 
Hyinc ignores the LCP and propagates through a vacuum, and then the 
beam of interest with Ex, Ey, Ez, Hx, and Hy. (Hz is identically zero.) 
It is a 1-dimensional problem with the beam propagating in the z- 
direction. 

9 - 9 ' 2 ' 2 ' 2 - 9 ' 2 ' 9 ' 2 - 9 ' 2 ' 9 ' 2 - 9 ' 5 ' 9 ' 2 - 2 ' 2 ' 9 ' 5 - 9 ' 2 ' 9 ' 2 - 9 ' 2 ' 9 ' 5 - 9 ' 2 - 9 ' 2 - 9 ' 2 - 2 ' 5 - 9 ' 2 - 9 ' 2 ' 2 ' 2 ' 9 ' 5 - 9 ' 2 - 2 ' 2 - 9 - 2 ' 2 ' 2 - 9 ' 2 - 9 ' 2 ' 9 ' 2 - 9 - 2 - 9 - 9 - 9 - 2 - 9 ' 2 ' 9 - 9 - 9 - 2 - 

ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

clear; 

%% Material parameters and anchoring conditions 
global h psibot thetabot psitop thetatop 

% Physical dimensions and the medium 
"6 

% The LCP layer has the physical width realh, which is chosen to be 10 
free space wavelengths of the incident beam. But, I nondimensionalize 
the system so that the width of the LCP layer is h=l. The anisotropic 
LCP layer goes from z=0 to z=l. The LCP is sandwiched between two 
isotropic glass layers of (nondimensional) thickness hglass. (However, 
for this code I wanted to ignore the reflections of the glass to focus 
on the effect of the LCP, and so I chose the glass to be optically the 
same as the air/vacuum, but to make the code work right, I still 
include the glass layers as one numerical cell thick.) Outside the 
glass layers are layers of air (really vacuum) of thickness hair. This 
is the physical system, but there are two more artificial layers used 
to help handle the difficulties of artificially truncating the 
computational domain without creating nonphysical reflections. The 
scattering layer is width hscat, and the perfectly matched layer (PML) 
is of width hPML. 

realh=6.330e-6; %physical gap between glass plates 

h=l; %rescaled gap width 

%hglass=h/4; %rescaled width of glass plates 
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hair=h/4; 
hscat=h/32; 
hPML=h/32; 


%rescaled width of the air layers 

%rescaled width of the scattering layer outside the PML 
%rescaled width of the PML 

epsilon_perp=l.5 A 2; %ordinary relative permittivity of the LCP 

epsilon_para=l.7 A 2 ; %extraordinary relative permittivity of the LCP 

delta_epsilon=epsilon_para-epsilon_perp; 

epsilon_air=l; %relative permittivity in air and scatter field 

epsilon_glass=l; %1.5 A 2; %relative permittivity in glass 

epsilon_perp 

% The incident beam takes the form 

% Exinc=EO*sin(wavenumber*(z-zstart)-omega*t)*(1-exp(-(t- 
tdelay) A 2/sigmadelay)) 

"6 

% The Gaussian decay term is there to ease the incident beam into the 
system. If it suddenly turns on at full force, it can introduce 
artificial reflections. I included the possibility of a tdelay, but I 
set it to zero. This decays to a simple sine wave. Hyinc is the same. 
Also, I have used the nondimensional timescale tO=realh/cO. 

E0=1; %strength of the incident field 

lambda=6.330e-7; %meters 

cO=l; %299792458; %speed of light in meters per second 
omega=2*pi*realh/lambda; %2*pi*c0/lambda; 

wavenumber=2*pi*realh/lambda; 
tdelay=0; 
sigmadeiay=0.5; 

EndTime=150*2*pi/omega % The final time 

% Description of the LCP orientation 
"6 

% This uses the Leslie-Ericksen (uniaxial) description of the 

orientation with a unit vector called the major director 

n= (cos (theta(z))*cos(psi(z)),cos(theta(z))*sin(psi(z)),sin(theta(z))) A T 

'O 

% The angles theta(z) and psi(z) are solutions to a system of ODEs that 
are re-solved at each successively finer grid, but the anchoring 
(boundary) conditions are our control variables in the system unixial 
anchoring conditions at the bottom of the LCP layer (z=0) 
psibot=2*pi/3; 

thetabot=pi/4; % theta=0 is tangential 

% unixial anchoring conditions at the top of the LCP layer (z=l) 
psitop=2*pi/3; 

thetatop=pi/4; % theta=0 is tangential 

ntheta=sqrt(epsilon_perp*epsilon_para/(epsilon_perp*(cos(thetabot)) A 2+e 
psilon_para*(sin(thetabot)) A 2)); % extraordinary index of refraction 
nperp=sqrt(epsilon_perp); % ordinary index of refraction 

kperp=omega*nperp/cO; % ordinary wave number 
ktheta=omega*ntheta/cO; % extraordinary wave number 
tic 

%% The coarsest grid 
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% Numerical parameters 

NLCP=100; % number of numerical cells in the LCP 

dz=h/NLCP; % width of one cell 

% The so-called CFL condition is a relationship between dt and dz that 
must 

% be satisfied by a grid in order for the solution to converge. In 
this 

% case, the CFL condition is cO*dt/dz <= 1. 
cdtoverdz=l; 

dt=dz*cdtoverdz; % time step 

Tfin=round(EndTime/dt)+1; % must have one more time step beyond 

EndTime to get H at Tfin 

timeperiod=round(2*pi/omega/dt); % the number time steps in one time 

period 


% The physical parts of the system (LCP, glass, and air layers) are 
called the total field. The other parts are part of the scattered 
field. In the total field, electromagnetic field variables E, D, and H 
store the sum of the scattered field and the incident beam, which we 
know at all locations and times; however, in the scattered field, these 
variables store only the scattered field. We do this because we do not 
want the source of the incident beam to be a part of the "real" system. 
By using a "hard source" (a point at which we specify what the 
electric field is), we create a point that makes an artificial 
reflection if the wave returns there. 

The PML uses two parameters to quantify how it absorbs the 
electromagnetic field without reflection 

sigmamax=2/dt; % in the PML, maximum value of sigma at the 

PEC boundary 

sigmam=3; % in the PML, sigmaz=sigmamax((kPMLbot+1- 

k)/kPMLbot) A sigmam 


Nglass = l; %round(hglass/dz) ; %number of numerical cells in the glass 
Nair=round(hair/dz) ; % number of numerical cells in the air 

between glass and scattering layer 

Nscat=round(hscat/dz); % number of numerical cells in the 

scattering later 

NPML=round(hPML/dz); % number of numerical cells in the 

perfectly matched layer 

kPMLbot=NPML; % sigma defined but equal to zero at 

k=kPMLbot 


kscatbot=kPMLbot+Nscat; 
total field; kscatbot-1/2 
kairbot=kscatbot+Nair; 
bottom plate 

kglassbot=kairbot+Nglass; 
plate/where the anchoring 
kglasstop=kglassbot+NLCP; 
plate/where the anchoring 
kairtop=kglasstop+Nglass; 
plate 

kscattop=kairtop+Nair; 
kscattop+1/2 is not 
kPMLtop=kscattop+Nscat; 
k=kPMLtop 


% k=kscatbot is considered to be in the 
is in the scattered field 

% k=kairbot is the bottom edge of the 

% k=kglassbot is the top edge of the bottom 
condition is applied 

% k=kglasstop is the bottom edge of the top 
condition is applied 

% k=kairtop is the top edge of the top 

% k=kscattop is in the total field; 

% sigma defined but equal to zero at 
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Ntot=kPMLtop+NPML; 
kincstart=kscatbot-2; 
ze=2:Ntot; 

in the D and E equations 
zh=2:Ntot + 1; 
in the H equations 


% Ntot cells in total 

% location of the incident field 

% k values for updating the gradients of H 


% k values for updating the gradients of E 


% Define electromagnetic field variables 


% This code is set up to store the electromagnetic variables for all 
values of z_k=(k-kglassbot-1)*dz but to only store the values of 
t_n=n*dt for the last time period (including both the first and 
lasttime of the period). We overwrite the values at each time step 
until we get to n=Tbreak 
Tbreak=Tfin-2*timeperiod-l; 

%at which time we switch to storing all the time steps. E and D are 
known on the space-time grid. At time zero, they are zero, but this is 
not stored. Ex and Dx are zero at k=l and k=Ntot+l but they are stored 
so that Ex(k,n) stores Ex at z=(k-kglassbot-1)*dz and t=n*dt, or Ex_k A n 
is stored as Ex(k+l,n). 

Dx=zeros(Ntot + 1,2*timeperiod+l) ; 

Dy=zeros(Ntot + 1,2*timeperiod+l) ; 

Dz=zeros(Ntot + 1,2*timeperiod+l) ; 

Ex=zeros(Ntot + 1,2*timeperiod+l) ; 

Ey=zeros(Ntot + 1,2*timeperiod+l) ; 

Ez=zeros(Ntot + 1,2*timeperiod+l) ; 

% H is known on the half-grid. Hx(k,n) is Hx at z=(k-kglassbot+1/2)*dz 
% and t= (n-1/2)*dt, or Hx_(k+1/2) A (n+1/2) is stored as Hx(k+l,n+l). 

% Hz is zero and not stored. We have to store one more time step so 
that we can average the half-grid times to get the on-grid values. 
Similarly, we average the half-grid spaces to get the on-grid values, 
but because we E and D store the zero values at the boundary, there is 
one more stored value of Ex than Hy. 

Hx=zeros(Ntot,2*timeperiod+2); 

Hy=zeros(Ntot,2*timeperiod+2) ; 

% The incident beam is stored as Exinc, Eyinc, Hxinc, and Hyinc 
although the incident beam we use has Eyinc=Hxinc=0. 

Exinc=zeros(Ntot + 1,2*timeperiod+l) ; 

Eyinc=zeros(Ntot + 1,2*timeperiod+l) ; 

Hxinc=zeros(Ntot,2*timeperiod+2 ) ; 

Hyinc=zeros(Ntot,2*timeperiod+2) ; 

% The "tilde" H vales interpolate the half-grid values of H onto the 
grid 

Hxtilde=zeros(Ntot-1,2*timeperiod+l); 

Hytilde=zeros(Ntot-1,2*timeperiod+l) ; 

Hxinctilde=zeros(Ntot-1,2*timeperiod+l) ; 

Hyinctilde=zeros(Ntot-1,2*timeperiod+l); 
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% The Poynting vector, the cross product of E and H, is known on the 
space-time grid. It is not stored at k=l or k=Ntot+l 
Poyntingx=zeros(Ntot-1,2*timeperiod+l) ; 

Poyntingy=zeros(Ntot-1,2*timeperiod+l) ; 

Poyntingz=zeros(Ntot-1,2*timeperiod+l) ; 

PoyntingMag=zeros(Ntot-1,2*timeperiod+l); %The magnitude of the 
Poynting vector 

% Compute the permittivity tensor epsilon 

~o 

% First, we must compute the LCP orientation by solving an ODE. It 

uses different z values used for computing the orientation 

z1=1inspace(0,h,NLCP+1); 

solinit=bvpinit(zl,@LEodeinit); 

sol=bvp4c(@LEodefun,@LEbcfun,solinit); 

u=deval(sol,zl); 

psiLE=u (1, : ); 

thetaLE=u (3, : ); 

nl=cos(thetaLE).*cos(psiLE); 
n2=cos(thetaLE).*sin(psiLE); 
n3=sin(thetaLE); 

% The permittivity tensor epsilon is symmetric and stored on the grid 
for the whole domain. In the LCP, we use the above orientation, but we 
also need to define it in the isotropic regions 

epsxx=zeros(Ntot+1,1); 
epsyy=zeros(Ntot+1,1); 
epszz=zeros(Ntot+1,1); 
epsxy=zeros(Ntot+1,1); 
epsxz=zeros(Ntot+1,1); 
epsyz=zeros(Ntot+1,1); 

epsxx(1:kPMLbot+1)=1; epsyy(1:kPMLbot+1)=1; epszz(1:kPMLbot+1)=1; 

epsxx(kPMLbot+2:kairbot)=epsilon_air; 

epsyy(kPMLbot+2:kairbot)=epsilon_air; 

epszz(kPMLbot+2:kairbot)=epsilon_air; 

epsxx(kairbot + 1:kglassbot)=epsilon_glass; 

epsyy(kairbot+1:kglassbot)=epsilon_glass; 

epszz (kairbot + 1:kglassbot)=epsilon_glass; 

epsxx(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*nl. A 2; 

epsyy(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n2. A 2; 

epszz(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n3. A 2; 

epsxy(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n2; 

epsxz(kglassbot + 1:kglasstop+1)=delta_epsilon*nl.*n3; 

epsyz(kglassbot+1:kglasstop+1)=delta_epsilon*n2.*n3; 

epsxx(kglasstop+2:kairtop+1)=epsilon_glass; 

epsyy(kglasstop+2:kairtop+1)=epsilon_glass; 

epszz(kglasstop+2:kairtop+1)=epsilon_glass; 

epsxx(kairtop+2:kPMLtop)=epsilon_air; 

epsyy(kairtop+2:kPMLtop)=epsilon_air; 

epszz(kairtop+2:kPMLtop)=epsilon_air; 
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epsxx(kPMLtop+1:Ntot)=1; epsyy(kPMLtop+1:Ntot)= 1; 
epszz(kPMLtop+1:Ntot)= 1; 

% For computation, we really need the inverse of epsilon, and so we 
will pre-compute that. 

epsdet=epsxx.*(epsyy.*epszz-epsyz. A 2)-epsxy.*(epsxy.*epszz- 

epsxz.*epsyz)+epsxz.*(epsxy.*epsyz-epsyy.*epsxz); 

epsinvxx=(epsyy.*epszz-epsyz. A 2) ./epsdet; 

epsinvxy=(epsxz.*epsyz-epsxy.*epszz)./epsdet; 

epsinvxz=(epsxy.*epsyz-epsyy.*epsxz)./epsdet; 

epsinvyy=(epsxx.*epszz-epsxz. A 2) ./epsdet; 

epsinvyz=(epsxy.*epsxz-epsxx.*epsyz)./epsdet; 

epsinvzz=(epsxx.*epsyy-epsxy. A 2)./epsdet; 


% Define coefficients used in our 
% This is where the PML absorbing 
C_E=ones(Ntottl, 1); 
update equations 

D_E=cdtoverdz*ones(Ntot + l, 1) ; % 

update equations 

C_H=ones(Ntot, 1) ; % 

update equations 

D_H=cdtoverdz*ones(Ntot, 1) ; % 

update equations 

for k=l:kPMLbot % 

PML 


update scheme. 

factors are used. 

factor of the D A n term in the D 

factor of the curl H term in the D 

factor of the H A n term in the H 

factor of the curl H term in the H 

fill in the effect of sigma in the 


sigmaz=sigmamax*((kPMLbot-k+1)/kPMLbot) A sigmam; 
C_E(k)=(2-dt*sigmaz)/(2tdt*sigmaz); 

D_E(k)=2*cdtoverdz/(2tdt*sigmaz) ; 

C_E(Ntot+2-k)=(2-dt*sigmaz)/(2tdt*sigmaz); 

D_E(Ntot+2-k)=2*cdtoverdz/(2+dt*sigmaz); 
sigmaz=sigmamax*((kPMLbot-k+0.5)/kPMLbot) A sigmam; 
C_H(k) = (2-dt*sigmaz)/(2+dt*sigmaz) ; 

D_H(k)=2*cdtoverdz/(2+dt*sigmaz) ; 

C_H(Ntottl-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_H(Ntottl-k)=2*cdtoverdz/(2tdt*sigmaz); 

end 


% The starting main loop, not storing every step 
Exinc(kincstarttl, 1)=-E0*sin(omega*dt)*(1-exp(-(dt- 
tdelay) A 2/sigmadelay)) ; 

Hx(kscatbot,1)=-D_H(kscatbot)*Eyinc(kscatbottl, 1) ; 

Hy(kscatbot,1)=D_H(kscatbot)*Exinc(kscatbottl,1); 

Hxinc(:,1)=D_H .* (Eyinc(zh,1)-Eyinc(zh-1,1)) ; 

Hyinc(:,1)=-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)) ; 
for n=2:Tbreak-l 

% Update dD/dt=curl(H) at t=(ntl)*dt 

Dx(ze,1)=C_E(ze) .*Dx(ze,1)-D_E(ze) .*(Hy(ze, 1)-Hy(ze-1,1)) ; 

Dy (ze, 1) =C_E (ze) . *Dy (ze, 1) tD_E (ze) . * (Hx (ze, 1) -Hx (ze-1, 1) ) ; 

% There is a correction at scattered field-total field boundary 

Dx(kscatbottl,1)=Dx(kscatbottl,1)tD_E(kscatbottl)*Hyinc(kscatbot,1); 

Dy(kscatbottl,1)=Dy(kscatbottl, 1) - 
D_E(kscatbottl)*Hxinc(kscatbot, 1) ; 


68 


Dx(kscattop+1,1)=Dx(kscattop+1,1)- 
D_E(kscattop+1)*Hyinc(kscattop+1,1); 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1, 1) ; 

% Get E from E=epsilon A -l*D 

Ex(ze,1)=epsinvxx(ze) .*Dx(ze, 1)+epsinvxy(ze) .*Dy(ze, 1)+epsinvxz(ze) .*Dz 
(ze,1) ; 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze,1)+epsinvyy(ze) .*Dy (ze,1)+epsinvyz(ze) .*Dz 
(ze,1) ; 

% Update the incident electric field 

Exinc(ze,1)=C_E(ze).*Exinc(ze,1)-D_E(ze).*(Hyinc(ze,1)-Hyinc(ze- 

i,D); 

Eyinc(ze,1)=C_E(ze).*Eyinc(ze,1)+D_E(ze).*(Hxinc(ze,1)-Hxinc(ze- 

i,D); 

Exinc(kincstart + 1,1)=-E0*sin(omega*n*dt)*(1-exp(-(n*dt- 
tdelay) A 2/sigmadelay)) ; 

% Update dH/dt=-curl(E) at t=(n+3/2)*dt 

Hx(:,1)=C_H.*Hx(:,1)+D_H.*(Ey(zh, 1)-Ey(zh-1,1)) ; 

Hy(:,1)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1,1)); 

% Correction at the scattered field-total field boundary 

Hx(kscatbot,1)=Hx(kscatbot,1)-D_H(kscatbot)*Eyinc(kscatbot+1,1); 
Hy(kscatbot,1)=Hy(kscatbot,1)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,1)=Hx(kscattop+1,1)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 

Hy (kscattop+1,1)=Hy(kscattop+1,1)- 
D_H(kscattop+1)*Exinc(kscattop+1, 1) ; 

% Update the incident magnetic field 

Hxinc(:,1)=C_H.*Hxinc(:,1)+D_H.*(Eyinc(zh,1)-Eyinc(zh-1,1)); 

Hyinc(:,1)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)) ; 

end 

% Now for the last period, we store each time step. Also, we need to 
start computing the Poynting vector, and so we need to compute Htilde, 
which is H interpolated onto the grid. The first time through is 
different than the following loop and the last time step. 

Dx(ze,1)=C_E(ze) .*Dx(ze,1)-D_E(ze) .*(Hy(ze, 1)-Hy(ze-1,1)) ; 

Dy(ze,1)=C_E(ze).*Dy(ze,1)+D_E(ze).*(Hx(ze,1)-Hx(ze-1,1)); 

Dx(kscatbot + 1, 1)=Dx(kscatbot + 1, 1)+D_E(kscatbot + 1)*Hyinc(kscatbot,1); 
Dy(kscatbot+1,1)=Dy(kscatbot+1,1)-D_E(kscatbot+1)*Hxinc(kscatbot,1); 
Dx (kscattop+1,1)=Dx(kscattop+1, 1)- 
D_E(kscattop+1)*Hyinc(kscattop+1,1); 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1); 
Ex(ze,1)=epsinvxx(ze).*Dx(ze,1)+epsinvxy(ze).*Dy(ze,1)+epsinvxz(ze).*Dz 
(ze,1) ; 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze,1)+epsinvyy(ze) .*Dy (ze,1)+epsinvyz(ze) .*Dz 
(ze,1) ; 

Ez(ze,1)=epsinvxz(ze).*Dx(ze,1)+epsinvyz(ze).*Dy(ze,1)+epsinvzz(ze). *Dz 
(ze,1); 

Exinc(ze,1)=C_E (ze) .*Exinc(ze,1)-D_E (ze) .*(Hyinc(ze,1)-Hyinc (ze-1,1)); 
Eyinc(ze,1)=C_E(ze) .*Eyinc(ze,1)+D_E(ze) .*(Hxinc(ze,1)-Hxinc (ze-1,1)); 
Exinc(kincstart+1,1)=-E0*sin(omega*Tbreak*dt)*(1-exp(-(Tbreak*dt- 
tdelay) A 2/sigmadelay)) ; 

Hx ( : , 2) =C_H. *Hx ( : , 1) +D_H. * (Ey (zh, 1) -Ey (zh-1,1) ) ; 
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Hy(:,2)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1,1)); 

Hx(kscatbot,2)=Hx(kscatbot,2)-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,2)=Hy(kscatbot,2)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,2)=Hx(kscattop+1,2)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 

Hy (kscattop+1,2)=Hy(kscattop+1,2)- 
D_H(kscattop+1)*Exinc(kscattop+1,1) ; 

Hxinc(:,2)=C_H.*Hxinc(:,1)+D_H.* (Eyinc(zh, 1)-Eyinc(zh-1, 1)) ; 

Hyinc(:,2)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh,1)-Exinc(zh-1, 1)) ; 

Hxinctilde(:,!) = (Hxinc(ze,1)+Hxinc(ze,2)+Hxinc(ze-1,1)+Hxinc (ze- 
1,2))/4; 

Hyinctilde(: , 1) = (Hyinc(ze, 1)+Hyinc(ze, 2 ) +Hyinc(ze-1,1)+Hyinc(ze- 
1,2))/4; 

Hxtilde(:,1)=(Hx(ze,1)+Hx(ze, 2 ) +Hx(ze-1,1)+Hx(ze-1,2))/4; 

Hytilde(:,1) = (Hy(ze,1)+Hy(ze, 2)+Hy(ze-1, 1)+Hy(ze-l,2))/4; 

Poyntingx(:,1)=-Ez(ze, 1) .*Hytilde(:,1); 

Poyntingy(:,l)=Ez(ze,l).*Hxtilde(:,1); 

Poyntingz(: , 1)=Ex(ze, 1) .*Hytilde(:,1)-Ey(ze,1) .*Hxtilde(:,1); 
PoyntingMag(:,1)=sqrt(Poyntingx(:,1).''2+Poyntingy(:,1).^2+Poyntingz(:,1 
). A 2); 

for n=l:2*timeperiod 

Dx(ze, n+1)=C_E(ze) .*Dx(ze, n)-D_E(ze) .*(Hy(ze,n+1)-Hy(ze-l,n+l)); 

Dy(ze,n+1)=C_E(ze).*Dy(ze,n)+D_E(ze).*(Hx(ze,n+1)-Hx(ze-1,n+1)); 

Dx(kscatbot+1,n+1)=Dx(kscatbot+1,n+1)+D_E(kscatbot+1)*Hyinc(kscatbot,n+ 

i); 

Dy (kscatbot + 1,n+1)=Dy(kscatbot + 1, n+1) - 
D_E(kscatbot + 1)*Hxinc(kscatbot, n+1); 

Dx(kscattop+1,n+1)=Dx(kscattop+1, n+1)- 
D_E(kscattop+1)*Hyinc(kscattop+1, n+1) ; 

Dy(kscattop+1,n+1)=Dy(kscattop+1,n+1)+D_E(kscattop+1)*Hxinc(kscattop+1, 
n+1) ; 

Ex(ze,n+1)=epsinvxx(ze).*Dx(ze,n+1)+epsinvxy(ze).*Dy(ze,n+1)+epsinvxz(z 
e).*Dz(ze,n+1); 

Ey(ze, n+1)=epsinvxy(ze) .*Dx(ze, n+1)+epsinvyy(ze) .*Dy(ze, n+1)+epsinvyz(z 
e).*Dz(ze,n+1); 

Ez(ze, n+1)=epsinvxz(ze) .*Dx(ze, n+1)+epsinvyz (ze) .*Dy(ze,n+1)+epsinvzz (z 
e).*Dz(ze,n+1); 

Exinc(ze,n+1)=C_E (ze) .*Exinc(ze,n)-D_E(ze) .* (Hyinc(ze,n+1)-Hyinc (ze- 
1,n+1)); 

Eyinc(ze, n+1)=C_E(ze) .*Eyinc(ze,n)+D_E(ze) .*(Hxinc(ze,n+1)-Hxinc (ze- 
1,n+1)); 

Exinc(kincstart+1,n+1)=-E0*sin(omega*(Tbreak+n)*dt)*(1-exp(- 
( (Tbreak+n)*dt-tdelay) A 2/sigmadelay)) ; 

Hx(:,n+2)=C_H.*Hx(:,n+1)+D_H.*(Ey(zh,n+1)-Ey(zh-1,n+1)); 

Hy(:,n+2)=C_H.*Hy(:,n+1)-D_H.*(Ex(zh,n+1)-Ex(zh-1, n+1)); 

Hx (kscatbot,n+2)=Hx(kscatbot, n+2)- 
D_H(kscatbot)*Eyinc(kscatbot + 1, n+1) ; 

Hy(kscatbot,n+2)=Hy(kscatbot,n+2)+D_H(kscatbot)*Exinc(kscatbot+1,n+1); 
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Hx(kscattop+1, n+2)=Hx(kscattop+1,n+2)+D_H(kscattop+1)*Eyinc(kscattop+1, 
n+1) ; 

Hy (kscattop+1,n+2)=Hy(kscattop+1, n+2) - 
D_H(kscattop+1)*Exinc(kscattop+1, n+1); 

Hxinc(:,n+2)=C_H.*Hxinc(:,n+1)+D_H.*(Eyinc(zh, n+1)-Eyinc(zh-1, n+1)) ; 
Hyinc(:,n+2)=C_H.*Hyinc(:,n+1)-D_H.*(Exinc(zh,n+1)-Exinc(zh-1,n+1)); 
Hxinctilde(:,n+1)=(Hxinc(ze,n+1)+Hxinc(ze,n+2)+Hxinc(ze- 
1,n+1)+Hxinc(ze-l,n+2))/4; 

Hyinctilde(: , n+1) = (Hyinc(ze,n+1)+Hyinc(ze,n+2)+Hyinc(ze- 
1,n+1)+Hyinc(ze-l,n+2))/4; 

Hxtilde(: , n+1) = (Hx(ze, n+1)+Hx(ze, n+2)+Hx(ze-1, n+1)+Hx(ze-1, n+2))/4; 
Hytilde(:,n+1) = (Hy(ze,n+1)+Hy(ze,n+2)+Hy(ze-1,n+1)+Hy(ze-1, n+2))/4 ; 
Poyntingx(:,n+1)=-Ez(ze,n+1).*Hytilde(:,n+1); 

Poyntingy(:,n+1)=Ez(ze,n+1).*Hxtilde(:,n+1); 

Poyntingz(:,n+1)=Ex(ze,n+1).*Hytilde(:,n+1)- 
Ey(ze,n+1).*Hxtilde(:,n+1); 

PoyntingMag(: , n+1)=sqrt(Poyntingx(: , n+1) . A 2+Poyntingy(:,n+1) . A 2+Poyntin 

gz(:,n+1) . A 2) ; 

end 

% Now we compute the intensity, which is the time-average of the 
magnitude of the Poynting vector over one period. We use the trapezoid 
rule to approximate the integration 
Intensity=PoyntingMag(:,1); 
for n=2:2*timeperiod 

Intensity=Intensity+2*PoyntingMag(:,n); 

end 

Intensity=Intensity+PoyntingMag(:,2*timeperiod+l); 
Intensity=Intensity*dt/EO A 2/(4*pi/omega); 

% We are now done with this grid, and so we store the values to compare 
later with the finer grids. 

kPMLbotl=kPMLbot; kscatbotl=kscatbot; kairbotl=kairbot; 
kglassbotl=kglassbot; 

kglasstopl=kglasstop; kairtopl=kairtop; kscattopl=kscattop; 

kPMLtopl=kPMLtop; 

timeperiodl=timeperiod; 

Ntotl=Ntot; dzl=dz; dtl=dt; 

Exl=Ex; Eyl=Ey; Ezl=Ez; 

Hxl=Hx; Hyl=Hy; 

Exincl=Exinc; Eyincl=Eyinc; 

Hxincl=Hxinc; Hyincl=Hyinc; 

Hxinctildel=Hxinctilde; Hyinctildel=Hyinctilde; 

Hxtildel=Hxtilde; Hytildel=Hytilde; 

Poyntingxl=Poyntingx; Poyntingyl=Poyntingy; Poyntingzl=Poyntingz; 
PoynitngMagl=PoyntingMag; Intensityl=Intensity; 

Tfinl=Tfin; Tbreakl=Tbreak; 
kincstartl=kincstart; 

zTotEl=dz*((1:Ntot+1)-kglassbot-1); % locations of E, D, and 

epsilon, z=0 is bottom of the LCP, z=l is top 

zTotHl=dz*((1:Ntot)-1/2-kglassbot); % locations where H is known; 

not contained in previous refinement 

toe 

tic 
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%% Finer grid 

% We now repeat the above loop for a grid spacing that is one-half the 
previous grid in both space and time, 
refinement=2; 


NLCP=NLCP*refinement; 

%number 

of 

numerical 

cells 

in 

the 

LCP 

dz=h/NLCP; 

Nglass=Nglass*refinement; 

%number 

of 

numerical 

cells 

in 

the 

glass 

Nair=Nair*refinement; 

%number 

of 

numerical 

cells 

in 

the 

air 

between glass and scattering 
Nscat=Nscat*refinement; 

layer 

%number 

of 

numerical 

cells 

in 

the 


scattering later 
NPML=NPML*refinement; 

%number 

of 

numerical 

cells 

in 

the 



perfectly matched layer 
cdtoverdz=l; 
dt=dz*cdtoverdz; 

Tfin=floor(EndTime/dt)+1; 
timeperiod=round(2*pi/omega/dt); 

sigmamax=2/dt; %10 A 17; %in the PML, sigmaz = sigmamax ( (kPML- 

1)/(NPML-1)) A sigmam 
sigmam=3; 

kPMLbot=NPML; % sigma defined but equal to zero at 

k=kPMLbot 


kscatbot=kPMLbot+Nscat; 
total field; kscatbot-1/2 
kairbot=kscatbot+Nair; 
bottom plate 

kglassbot=kairbot+Nglass; 
plate/where the anchoring 
kglasstop=kglassbot+NLCP; 
plate/where the anchoring 
kairtop=kglasstop+Nglass; 
plate 

kscattop=kairtop+Nair; 
kscattop+1/2 is not 
kPMLtop=kscattop+Nscat; 
k=kPMLtop 
Ntot=kPMLtop+NPML; 
ze=2:Ntot; 

in the D and E equations 
zh=2:Ntot+1; 

z2 = linspace(0,h,NLCP + 1) ; 
kincstart=refinement*(kinc 


% k=kscatbot is considered to be in the 
is in the scattered field 

% k=kairbot is the bottom edge of the 

% k=kglassbot is the top edge of the bottom 
condition is applied 

% k=kglasstop is the bottom edge of the top 
condition is applied 

% k=kairtop is the top edge of the top 

% k=kscattop is in the total field; 

% sigma defined but equal to zero at 

% Ntot cells in total% ze=2:Ntot; 

% k values for updating the gradients of H 


%z values used for computing the orientation 
startl-kglassbotl)+kglassbot; 


%% compute the permittivity 
solinit=bvpinit(z2,@LEodeinit) ; 
sol=bvp4c(@LEodefun,@LEbcfun, solinit) ; 
u=deval(sol, z2) ; 
psiLE=u (1, : ); 
thetaLE=u (3, : ); 


nl=cos(thetaLE).*cos(psiLE); 
n2=cos(thetaLE).*sin(psiLE); 
n3=sin(thetaLE); 


epsxx=zeros(Ntot + 1,1) ; 
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e P s yy =zeros (Ntot + l, 1) ; 
epszz=zeros(Ntot+l,1); 
epsxy=zeros(Ntot + l, 1) ; 
epsxz=zeros(Ntot + l, 1) ; 
epsyz=zeros(Ntot+l,1); 


epsxx(1:kPMLbot + 1)=1; epsyy(1:kPMLbot + 1)=1; epszz(1:kPMLbot + 1)=1; 

epsxx(kPMLbot+2:kairbot)=epsilon_air; 

epsyy(kPMLbot+2:kairbot)=epsilon_air; 

epszz(kPMLbot+2:kairbot)=epsilon_air; 

epsxx(kairbot + 1:kglassbot)=epsilon_glass; 

epsyy(kairbot + 1:kglassbot)=epsilon_glass; 

epszz(kairbot + 1:kglassbot)=epsilon_glass; 

epsxx(kglassbot + 1:kglasstop+1)=epsilon_perp+delta_epsilon*nl. A 2; 

epsyy(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n2. A 2; 

epszz(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n3. A 2; 

epsxy(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n2; 

epsxz(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n3; 

epsyz(kglassbot+1:kglasstop+1)=delta_epsilon*n2.*n3; 

epsxx(kglasstop+2:kairtop+1)=epsilon_glass; 

epsyy(kglasstop+2:kairtop+1)=epsilon_glass; 

epszz(kglasstop+2:kairtop+1)=epsilon_glass; 

epsxx(kairtop+2:kPMLtop)=epsilon_air; 

epsyy(kairtop+2:kPMLtop)=epsilon_air; 

epszz(kairtop+2:kPMLtop)=epsilon_air; 

epsxx(kPMLtop+1:Ntot)=1; epsyy(kPMLtop+1:Ntot)=1; 

epszz(kPMLtop+1:Ntot)=1; 


epsdet=epsxx.*(epsyy.*epszz-epsyz. A 2)-epsxy.*(epsxy.*epszz- 

epsxz.*epsyz)+epsxz.*(epsxy.*epsyz-epsyy.*epsxz) ; 

epsinvxx=(epsyy.*epszz-epsyz. A 2)./epsdet; 

epsinvxy=(epsxz.*epsyz-epsxy.*epszz)./epsdet; 

epsinvxz=(epsxy.*epsyz-epsyy.*epsxz)./epsdet; 

epsinvyy=(epsxx.*epszz-epsxz. A 2)./epsdet; 

epsinvyz=(epsxy.*epsxz-epsxx.*epsyz)./epsdet; 

epsinvzz=(epsxx.*epsyy-epsxy. A 2)./epsdet; 


% Define coefficients 

C_E=ones(Ntot + l, 1) ; 

update equations 

D_E=cdtoverdz*ones(Ntot + l, 1) ; 

update equations 

C_H=ones(Ntot,1); 

update equations 

D_H=cdtoverdz*ones(Ntot, 1) ; 

update equations 

for k=l:kPMLbot 

PML 


% factor of the D A n term in the D 

% factor of the curl H term in the D 

% factor of the H A n term in the H 

% factor of the curl H term in the H 

% fill in the effect of sigma in the 


sigmaz=sigmamax*((kPMLbot-k+1)/kPMLbot) A sigmam; 
C_E(k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(k)=2*cdtoverdz/(2+dt*sigmaz); 

C_E(Ntot+2-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(Ntot+2-k)=2*cdtoverdz/(2+dt*sigmaz); 
sigmaz=sigmamax*((kPMLbot-k+0.5)/kPMLbot) A sigmam; 
C_H(k)=(2-dt*sigmaz)/(2+dt*sigmaz); 
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D_H(k)=2*cdtoverdz/(2+dt*sigmaz) ; 

C_H(Ntot+l-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 
D_H(Ntot+l-k)=2*cdtoverdz/(2+dt*sigmaz); 

end 


% Define electromagnetic field variable 
Dx=zeros(Ntot + 1,2*timeperiod+l) ; % 

time grid 

Dy=zeros(Ntot + 1,2*timeperiod+l) ; % 

initially (n=0), but this is not stored. 
Dz=zeros(Ntot+1,2*timeperiod+l); % 

k=Ntot+l but they are stored so that 
Ex=zeros(Ntot+1,2*timeperiod+l); % 

z=(k-kglassbot-1)*dz, 

Ey=zeros(Ntot + 1,2*timeperiod+l) ; % 

Ex(k+1,n) 

Ez=zeros(Ntot + 1,2*timeperiod+l) ; % 

Hx=zeros(Ntot,2*timeperiod+2) ; % 

z=(k-kglassbot+1/2)*dz 

Hy=zeros(Ntot,2*timeperiod+2) ; % 

Hx(k+1,n+1) 

Exinc=zeros(Ntot + 1,2*timeperiod+l) ; % 

can be nonzero, dependending on epsilon 
Eyinc=zeros(Ntot + 1,2*timeperiod+l) ; % 

Hxinc=zeros(Ntot,2*timeperiod+2) ; % 

are the incident field 
Hyinc=zeros(Ntot,2*timeperiod+2) ; 
Hxinctilde=zeros(Ntot-1,2*timeperiod+l) ; 
Hxinc onto the grid 

Hyinctilde=zeros(Ntot-1,2*timeperiod+l) ; 
Hxtilde=zeros(Ntot-1,2*timeperiod+l); 
onto the grid 

Hytilde=zeros(Ntot-1,2*timeperiod+l); 
Poyntingx=zeros(Ntot-1,2*timeperiod+l) ; 
on the space-time grid, and so H must be 
Poyntingy=zeros(Ntot-1,2*timeperiod+l) ; 
space. Thus it is not defined at k=l or 
Poyntingz=zeros(Ntot-1,2*timeperiod+l) ; 
PoyntingMag=zeros(Ntot-1,2*timeperiod+l) 
Tbreak=Tfin-2*timeperiod-l; 


D and E are known on the space- 

D and E fields are zero 

Dx and Ex are zero at k=l and 

Dx(k,n) is Dx at t=n*dt and 

or Ex_k A n is stored as 

H is known on the half-grid. 
Hx(k,n) is Hx at t=(n-1/2)*dt, 

Hx_(k+1/2) A (n+1/2) is stored as 

Hz is zero and not stored. Ez 

Exinc, Eyinc, Hxinc, and Hyinc 

% Hxinctilde interpolates 

% Hxtilde interpolates Hx 

% The Poynting vector is known 

% interpolated in time and in 

k=Ntot+l 

% or at t=Tfin 


% The starting main loop, not storing every step 
Exinc(kincstart + 1,1)=-E0*sin(omega*dt)*(1-exp(-(dt- 
tdelay) A 2/sigmadelay)) ; 

Hx(kscatbot,1)=-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,1)=D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hxinc(:,1)=D_H .* (Eyinc(zh,1)-Eyinc(zh-1, 1)) ; 

Hyinc(:,1)=-D_H.*(Exinc(zh,1)-Exinc(zh-1,1) ) ; 
for n=2:Tbreak-l 

Dx (ze, 1) =C_E (ze) . *Dx (ze, 1) -D_E (ze) . * (Hy (ze, 1) -Hy (ze-1, 1) ) ; 

Dy(ze,1)=C_E(ze).*Dy(ze,1)+D_E(ze).*(Hx(ze,1)-Hx(ze-1,1)); 

Dx(kscatbot+1,1)=Dx(kscatbot+1,1)+D_E(kscatbot+1)*Hyinc(kscatbot,1); 

Dy(kscatbot + 1,1)=Dy(kscatbot + 1,1)- 
D_E(kscatbot + 1)*Hxinc(kscatbot, 1) ; 
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Dx(kscattop+1,1)=Dx(kscattop+1,1)- 
D_E(kscattop+1)*Hyinc(kscattop+1,1); 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1, 1) ; 

Ex(ze,1)=epsinvxx(ze) .*Dx(ze,1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz (ze) .*Dz 
(ze,1); 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze, 1)+epsinvyy (ze) .*Dy(ze, 1)+epsinvyz (ze) . *Dz 
(ze,1); 

Exinc(ze,1)=C_E (ze) .*Exinc(ze,1)-D_E(ze) .*(Hyinc(ze,1)-Hyinc (ze- 

i,D); 

Eyinc(ze,1)=C_E(ze).*Eyinc(ze,1)+D_E(ze).*(Hxinc(ze,1)-Hxinc(ze- 

i,D); 

Exinc(kincstart + 1,1)=-EO*sin(omega*n*dt)*(1-exp(-(n*dt- 
tdelay)^2/sigmadelay)); 

Hx(: , 1)=C_H.*Hx(:,1)+D_H.*(Ey(zh,1)-Ey (zh-1,1)); 

Hy(:,1)=C_H.*Hy(:,1)-D_H.*(Ex(zh, 1)-Ex(zh-1,1)) ; 

Hx(kscatbot,1)=Hx(kscatbot,1)-D_H(kscatbot)*Eyinc(kscatbot+1,1); 

Hy(kscatbot,1)=Hy(kscatbot,1)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,1)=Hx(kscattop+1,1)+D_H(kscattop+1)*Eyinc(kscattop+1, 1) ; 

Hy(kscattop+1,1)=Hy(kscattop+1,1)- 
D_H(kscattop+1)*Exinc(kscattop+1, 1) ; 

Hxinc(:,1)=C_H.*Hxinc(:,1)+D_H.*(Eyinc(zh,1)-Eyinc(zh-1,1)); 

Hyinc(:,1)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)) ; 

end 

Dx(ze,1)=C_E(ze) .*Dx(ze,1)-D_E(ze) .*(Hy(ze, 1)-Hy(ze-1,1)) ; 

Dy(ze,1)=C_E(ze).*Dy(ze,1)+D_E(ze).*(Hx(ze,1)-Hx(ze-1,1)); 

Dx(kscatbot+1,1)=Dx(kscatbot+1,1)+D_E(kscatbot+1)*Hyinc(kscatbot,1); 
Dy(kscatbot+1,1)=Dy(kscatbot+1,1)-D_E(kscatbot+1)*Hxinc(kscatbot,1); 
Dx (kscattop+1,1)=Dx(kscattop+1,1)- 
D_E(kscattop+1)*Hyinc(kscattop+1, 1) ; 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1); 
Ex(ze,1)=epsinvxx(ze) .*Dx(ze,1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz (ze) .*Dz 
(ze,1); 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze, 1)+epsinvyy (ze) .*Dy(ze, 1)+epsinvyz (ze) . *Dz 
(ze,1) ; 

Ez(ze, 1)=epsinvxz(ze) .*Dx(ze,1)+epsinvyz(ze) .*Dy (ze,1)+epsinvzz(ze) . *Dz 
(ze,1) ; 

Exinc(ze,1)=C_E(ze) .*Exinc(ze,1)-D_E(ze) .*(Hyinc(ze, 1)-Hyinc(ze-1, 1)) ; 
Eyinc(ze, 1)=C_E(ze) .*Eyinc(ze,1)+D_E (ze) . *(Hxinc(ze,1)-Hxinc(ze-1,1)); 
Exinc(kincstart + 1,1)=-E0*sin(omega*Tbreak*dt)*(1-exp(-(Tbreak*dt- 
tdelay) A 2/sigmadelay)); 

Hx ( : ,2)=C_H.*Hx ( : , 1) +D_H. * (Ey (zh, 1) -Ey (zh-1, 1) ) ; 

Hy(:,2)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1,1)); 

Hx(kscatbot,2)=Hx(kscatbot,2)-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,2)=Hy(kscatbot,2)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,2)=Hx(kscattop+1,2)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 

Hy (kscattop+1,2)=Hy(kscattop+1,2)- 
D_H(kscattop+1)*Exinc(kscattop+1, 1) ; 

Hxinc(:,2)=C_H.*Hxinc(:,1)+D_H.*(Eyinc(zh,1)-Eyinc(zh-1,1)); 

Hyinc(:,2)=C_H.*Hyinc(:,1)-D_H .* (Exinc(zh,1)-Exinc(zh-1,1)); 
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Hxinctilde(:,1)=(Hxinc(ze,1)+Hxinc(ze,2)+Hxinc(ze-1,1)+Hxinc(ze- 
1,2))/4; 

Hyinctilde(:,1) = (Hyinc(ze,1)+Hyinc(ze,2)+Hyinc(ze-1,1)+Hyinc (ze- 
1,2))/4; 

Hxtilde(:,1) = (Hx(ze,1)+Hx(ze,2)+Hx(ze-1,1)+Hx (ze-l,2))/4; 

Hytilde(:,l)=(Hy(ze,l)+Hy(ze,2)+Hy(ze-1,1)+Hy(ze-1,2))/4; 

Poyntingx(:,1)=-Ez(ze,1).*Hytilde(:,1); 

Poyntingy(:,l)=Ez(ze,l).*Hxtilde(:,1); 

Poyntingz(:,1)=Ex(ze,1).*Hytilde(:,1)-Ey(ze,1).*Hxtilde(:,1); 
PoyntingMag(: , 1)=sqrt(Poyntingx(: , 1) . A 2+Poyntingy(: , 1) . A 2+Poyntingz(:,1 
) . A 2) ; 

for n=l:2*timeperiod 

Dx(ze,n+1)=C_E(ze) .*Dx(ze,n)-D_E(ze) .*(Hy(ze, n+1)-Hy(ze-l,n+l)) ; 

Dy(ze,n+1)=C_E(ze).*Dy(ze,n)+D_E(ze). *(Hx(ze,n+1)-Hx(ze-1,n+1)); 

Dx(kscatbot + 1,n+1)=Dx(kscatbot + 1,n+1)+D_E(kscatbot + 1)*Hyinc(kscatbot, n+ 

i); 

Dy(kscatbot + 1,n+1)=Dy(kscatbot + 1, n+1)- 
D_E(kscatbot + 1)*Hxinc(kscatbot, n+1) ; 

Dx(kscattop+1,n+1)=Dx(kscattop+1, n+1)- 
D_E(kscattop+1)*Hyinc(kscattop+1, n+1) ; 

Dy(kscattop+1,n+1)=Dy(kscattop+1,n+1)+D_E(kscattop+1)*Hxinc(kscattop+1, 
n+1) ; 

Ex(ze, n+1)=epsinvxx(ze) .*Dx(ze,n+1)+epsinvxy(ze) .*Dy(ze,n+1)+epsinvxz(z 
e).*Dz(ze,n+1); 

Ey(ze, n+1)=epsinvxy(ze) .*Dx(ze, n+1)+epsinvyy(ze) .*Dy(ze, n+1)+epsinvyz(z 
e).*Dz(ze,n+1); 

Ez(ze,n+1)=epsinvxz(ze) .*Dx(ze,n+1)+epsinvyz(ze) .*Dy(ze,n+1)+epsinvzz (z 
e).*Dz(ze,n+1); 

Exinc(ze,n+1)=C_E (ze) .*Exinc(ze,n)-D_E(ze) . *(Hyinc(ze,n+1)-Hyinc (ze- 
1,n+1)); 

Eyinc(ze,n+1)=C_E(ze).*Eyinc(ze,n)+D_E(ze). *(Hxinc(ze,n+1)-Hxinc(ze- 
1,n+1)); 

Exinc(kincstart+1,n+1)=-E0*sin(omega*(Tbreak+n)*dt)*(1-exp(- 
( (Tbreak+n)*dt-tdelay) A 2/sigmadelay)) ; 

Hx(:,n+2)=C_H.*Hx(:,n+1)+D_H.*(Ey(zh,n+1)-Ey(zh-1,n+1)); 

Hy(:,n+2)=C_H.*Hy(:,n+1)-D_H.*(Ex(zh,n+1)-Ex(zh-1, n+1)); 

Hx (kscatbot,n+2)=Hx(kscatbot, n+2)- 
D_H(kscatbot)*Eyinc(kscatbot + 1, n+1) ; 

Hy(kscatbot,n+2)=Hy(kscatbot,n+2)+D_H(kscatbot)*Exinc(kscatbot+1,n+1); 

Hx(kscattop+1,n+2)=Hx(kscattop+1,n+2)+D_H(kscattop+1)*Eyinc(kscattop+1, 
n+1) ; 

Hy(kscattop+1,n+2)=Hy(kscattop+1, n+2)- 
D_H(kscattop+1)*Exinc(kscattop+1, n+1) ; 

Hxinc(:,n+2)=C_H.*Hxinc(:,n+1)+D_H.*(Eyinc(zh,n+1)-Eyinc(zh-1, n+1)) ; 
Hyinc(:,n+2)=C_H.*Hyinc(:,n+1)-D_H.*(Exinc(zh,n+1)-Exinc(zh-1, n+1)) ; 
Hxinctilde(: , n+1) = (Hxinc(ze,n+1)+Hxinc(ze,n+2)+Hxinc (ze- 
1,n+1)+Hxinc(ze-1,n+2))/4; 
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Hyinctilde(:,n+1)=(Hyinc(ze,n+1)+Hyinc(ze,n+2)+Hyinc(ze- 
1,n+1)+Hyinc(ze-l,n+2))/4; 

Hxtilde(:,n+1) = (Hx(ze,n+1)+Hx(ze,n+2)+Hx(ze-1,n+1)+Hx(ze-1, n+2))/4; 
Hytilde(: , n+1) = (Hy(ze,n+1)+Hy(ze,n+2)+Hy(ze-1,n+1)+Hy (ze-1,n+2))/4; 
Poyntingx(:,n+1)=-Ez(ze,n+1).*Hytilde(:,n+1); 

Poyntingy(:,n+1)=Ez(ze,n+1).*Hxtilde(:,n+1); 

Poyntingz(:,n+1)=Ex(ze,n+1).*Hytilde(:,n+1)- 
Ey(ze,n+1).*Hxtilde(:,n+1); 

PoyntingMag(:,n+1)=sqrt(Poyntingx(:,n+1). A 2+Poyntingy(:,n+1). A 2+Poyntin 

gz(:,n+1) . A 2) ; 

end 

Intensity=PoyntingMag(: , 1) ; 
for n=2:2*timeperiod 

Intensity=Intensity+2*PoyntingMag(: , n) ; 

end 

Intensity=Intensity+PoyntingMag(:,2*timeperiod+l); 
Intensity=Intensity*dt/EO A 2/(4*pi/omega); 


kPMLbot2=kPMLbot; kscatbot2=kscatbot; kairbot2=kairbot; 
kglassbot2=kglassbot; 

kglasstop2=kglasstop; kairtop2=kairtop; kscattop2=kscattop; 
kPMLtop2=kPMLtop; 

Ntot2=Ntot; dz2=dz; dt2=dt; 

Ex2=Ex; Ey2=Ey; Ez2=Ez; 

Hx2=Hx; Hy2=Hy; 

Exinc2=Exinc; Eyinc2=Eyinc; 

Hxinc2=Hxinc; Hyinc2=Hyinc; 

Hxinctilde2=Hxinctilde; Hyinctilde2=Hyinctilde; 

Hxtilde2=Hxtilde; Hytilde2=Hytilde; 

Poyntingx2=Poyntingx; Poyntingy2=Poyntingy; Poyntingz2=Poyntingz; 

PoynitngMag2=PoyntingMag; Intensity2=Intensity; 

timeperiod2=timeperiod; Tfin2=Tfin; Tbreak2=Tbreak; 

kincstart2=kincstart; 

zTotE2=dz*((1:Ntot+1)-kglassbot-1); 

zTotH2=dz*((1:Ntot)-1/2-kglassbot); 

toe 

tic 

%% Finer grid 
refinement=2; 

NLCP=NLCP*refinement %number of numerical cells in the LCP 


in the glass 
in the air 


dz=h/NLCP; 

Nglass=Nglass*refinement; %number of numerical cells 

Nair=Nair*refinement; %number of numerical cells 

between glass and scattering layer 

Nscat=Nscat*refinement; %number of numerical cells in the 

scattering later 

NPML=NPML*refinement; %number of numerical cells in the 

perfectly matched layer 
cdtoverdz=l; 
dt=dz*cdtoverdz; 

Tfin=floor(EndTime/dt)+1; 
timeperiod=round(2*pi/omega/dt); 

sigmamax=2/dt; %in the PML, sigmaz=sigmamax((kPML-1)/(NPML- 

1)) A sigmam 
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sigmam=3; 

kPMLbot=NPML; % sigma defined but equal to zero at 

k=kPMLbot 


kscatbot=kPMLbot+Nscat; 
total field; kscatbot-1/2 
kairbot=kscatbot+Nair; 
bottom plate 

kglassbot=kairbot+Nglass; 


% k=kscatbot is considered to be in the 
is in the scattered field 

% k=kairbot is the bottom edge of the 

% k=kglassbot is the top edge of the bottom 


plate/where the anchoring condition is applied 


k=kglasstop is the bottom edge of the top 
ition is applied 

k=kairtop is the top edge of the top 

k=kscattop is in the total field; 

sigma defined but equal to zero at 

Ntot cells in total% ze=2:Ntot; 
h k values for updating the gradients of H 


kglasstop=kglassbot+NLCP; 

plate/where the anchoring condition is applied 
kairtop=kglasstop+Nglass; 
plate 

kscattop=kairtop+Nair; 
kscattop+1/2 is not 
kPMLtop=kscattop+Nscat; 
k=kPMLtop 
Ntot=kPMLtop+NPML; 
ze=2:Ntot; 
in the D and E equations 
zh=2:Ntot+1; 

z3=linspace(0,h,NLCP+1); %z values used for computing the orientation 
kincstart=refinement*(kincstart2-kglassbot2)+kglassbot; 

%% compute the permittivity 
solinit=bvpinit(z3,@LEodeinit); 
sol=bvp4c(SLEodefun,@LEbcfun,solinit); 
u=deval(sol, z3) ; 
psiLE=u (1, : ); 
thetaLE=u (3, : ); 


nl=cos(thetaLE).*cos(psiLE); 
n2=cos(thetaLE).*sin(psiLE); 
n3=sin(thetaLE); 


epsxx=zeros(Ntot + 1, 1) ; 
epsyy=zeros(Ntot + 1,1) ; 
epszz=zeros(Ntot+1,1); 
epsxy=zeros(Ntot + 1,1) ; 
epsxz = zeros(Ntot + 1,1) ; 
epsyz=zeros(Ntot+1,1); 

epsxx(1:kPMLbot+1)=1; epsyy(1:kPMLbot+1)=1; epszz(1:kPMLbot+1)=1; 

epsxx(kPMLbot+2:kairbot)=epsilon_air; 

epsyy(kPMLbot+2:kairbot)=epsilon_air; 

epszz(kPMLbot+2:kairbot)=epsilon_air; 

epsxx(kairbot+1:kglassbot)=epsilon_glass; 

epsyy(kairbot+1:kglassbot)=epsilon_glass; 

epszz (kairbot + 1:kglassbot)=epsilon_glass; 

epsxx(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*nl. A 2; 
epsyy(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n2. A 2; 
epszz(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n3. A 2; 
epsxy(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n2; 
epsxz(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n3; 
epsyz(kglassbot+1:kglasstop+1)=delta_epsilon*n2.*n3; 
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epsxx(kglasstop+2:kairtop+1)=epsilon_glass; 
epsyy(kglasstop+2:kairtop+1)=epsilon_glass; 
epszz(kglasstop+2:kairtop+1)=epsilon_glass; 
epsxx(kairtop+2:kPMLtop)=epsilon_air; 
epsyy(kairtop+2:kPMLtop)=epsilon_air; 
epszz(kairtop+2:kPMLtop)=epsilon_air; 
epsxx(kPMLtop+1:Ntot)=1; epsyy(kPMLtop+1:Ntot)=1; 
epszz(kPMLtop+1:Ntot)=1; 

epsdet=epsxx.*(epsyy.*epszz-epsyz. A 2)-epsxy.*(epsxy.*epszz- 

epsxz.*epsyz)+epsxz.*(epsxy.*epsyz-epsyy.*epsxz); 

epsinvxx=(epsyy.*epszz-epsyz."2)./epsdet; 

epsinvxy=(epsxz.*epsyz-epsxy.*epszz) ./epsdet; 

epsinvxz=(epsxy.*epsyz-epsyy.*epsxz)./epsdet; 

epsinvyy=(epsxx.*epszz-epsxz. A 2) ./epsdet; 

epsinvyz=(epsxy.*epsxz-epsxx.*epsyz)./epsdet; 

epsinvzz=(epsxx.*epsyy-epsxy. A 2)./epsdet; 

%% Define coefficients 


C_E=ones(Ntot + 1, 1); 

update equations 

D_E=cdtoverdz*ones(Ntot + 1, 1) ; 

update equations 

C_H=ones(Ntot,1); 

update equations 

D_H=cdtoverdz*ones(Ntot, 1); 

update equations 

for k=l:kPMLbot 

PML 


factor of the D A n term in the D 


factor of the curl H term in the D 


factor of the H A n term in the H 


factor of the curl H term in the H 


fill in the effect of sigma in the 


sigmaz=sigmamax*((kPMLbot-k+1)/kPMLbot) A sigmam; 
C_E(k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(k)=2*cdtoverdz/(2+dt*sigmaz) ; 

C_E(Ntot+2-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(Ntot+2-k)=2*cdtoverdz/(2+dt*sigmaz); 
sigmaz=sigmamax*((kPMLbot-k+0.5)/kPMLbot) A sigmam; 
C_H(k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_H(k)=2*cdtoverdz/(2+dt*sigmaz) ; 

C_H(Ntot+l-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_H(Ntot+l-k)=2*cdtoverdz/(2+dt*sigmaz); 


end 

%% Define electromagnetic field variables 


Dx=zeros(Ntot+1,2*timeperiod+l); % 

time grid 

Dy=zeros(Ntot + 1,2*timeperiod+l) ; % 

initially (n=0), but this is not stored 
Dz=zeros(Ntot+1,2*timeperiod+l); % 

k=Ntot+l but they are stored so that 
Ex=zeros(Ntot+1,2*timeperiod+l); % 

z=(k-kglassbot-1 )*dz, 

Ey=zeros(Ntot + 1,2*timeperiod+l) ; % 

Ex (k+1, n) 

Ez=zeros(Ntot + 1,2*timeperiod+l) ; % 

Hx=zeros(Ntot,2*timeperiod+2 ) ; % 

z=(k-kglassbot+1/2)*dz 

Hy=zeros(Ntot,2*timeperiod+2 ) ; % 

Hx(k+1 , n+1) 


D and E are known on the space- 

D and E fields are zero 

Dx and Ex are zero at k=l and 

Dx(k,n) is Dx at t=n*dt and 

or Ex_k A n is stored as 

H is known on the half-grid. 
Hx(k,n) is Hx at t=(n-1/2)*dt, 

Hx_(k+1/2) A (n+1/2) is stored as 
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Exinc=zeros(Ntot+1,2*timeperiod+l); % Hz is zero and not stored. Ez 

can be nonzero, dependending on epsilon 
Eyinc=zeros(Ntot+1,2*timeperiod+l); % 

Hxinc=zeros(Ntot,2*timeperiod+2); % Exinc, Eyinc, Hxinc, and Hyinc 

are the incident field 
Hyinc=zeros(Ntot,2*timeperiod+2); 

Hxinctilde=zeros(Ntot-1,2*timeperiod+l); % Hxinctilde interpolates 

Hxinc onto the grid 

Hyinctilde=zeros(Ntot-1,2*timeperiod+l); 

Hxtilde=zeros(Ntot-1,2*timeperiod+l); % Hxtilde interpolates Hx 

onto the grid 

Hytilde=zeros(Ntot-1,2*timeperiod+l); 

Poyntingx=zeros(Ntot-1,2*timeperiod+l) ; % The Poynting vector is known 

on the space-time grid, and so H must be 

Poyntingy=zeros(Ntot-1,2*timeperiod+l); % interpolated in time and in 

space. Thus it is not defined at k=l or k=Ntot+l 
Poyntingz=zeros(Ntot-1,2*timeperiod+l); % or at t=Tfin 

PoyntingMag=zeros(Ntot-1,2*timeperiod+l); 

Tbreak=Tfin-2*timeperiod-l; 

% The starting main loop, not storing every step 
Exinc(kincstart + 1,1)=-E0*sin(omega*dt)*(1-exp(-(dt- 
tdelay)^2/sigmadelay)) ; 

Hx(kscatbot,1)=-D_H(kscatbot)*Eyinc(kscatbot + 1,1) ; 

Hy(kscatbot,1)=D_H(kscatbot)*Exinc(kscatbot + 1,1) ; 

Hxinc(:,1)=D_H .* (Eyinc(zh,1)-Eyinc(zh-1, 1)) ; 

Hyinc(:,1)=-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)) ; 
for n=2:Tbreak-l 

Dx(ze,1)=C_E(ze).*Dx(ze,1)-D_E(ze).*(Hy(ze,1)-Hy(ze-1,1)); 

Dy (ze, 1) =C_E (ze) . *Dy (ze, 1) +D_E (ze) . * (Hx (ze, 1) -Hx (ze-1, 1) ) ; 

Dx(kscatbot + 1,1)=Dx(kscatbot + 1,1)+D_E(kscatbot + 1)*Hyinc(kscatbot, 1) ; 
Dy(kscatbot + 1,1)=Dy(kscatbot + 1,1)-D_E(kscatbot + 1)*Hxinc(kscatbot, 1) ; 
Dx(kscattop+1,1)=Dx(kscattop+1, 1) - 
D_E(kscattop+1)*Hyinc(kscattop+1,1); 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1) ; 

Ex(ze,1)=epsinvxx(ze) .*Dx(ze,1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz (ze) .*Dz 
(ze,1); 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze, 1)+epsinvyy (ze) .*Dy(ze, 1)+epsinvyz (ze) . *Dz 
(ze,1); 

Exinc(ze,1)=C_E (ze) .*Exinc(ze,1)-D_E(ze) .*(Hyinc(ze,1)-Hyinc (ze- 

i,D); 

Eyinc(ze,1)=C_E(ze) .*Eyinc(ze, 1)+D_E(ze) .*(Hxinc(ze,1)-Hxinc (ze- 

i,D); 

Exinc(kincstart+1,1)=-E0*sin(omega*n*dt)*(1-exp(-(n*dt- 
tdelay)^2/sigmadelay) ) ; 

Hx ( : , 1) =C_H. *Hx ( : , 1) +D_H. * (Ey (zh, 1) -Ey (zh-1,1) ) ; 

Hy(:,1)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1,1)); 

Hx(kscatbot,1)=Hx(kscatbot,1)-D_H(kscatbot)*Eyinc(kscatbot+1,1); 

Hy(kscatbot,1)=Hy(kscatbot,1)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,1)=Hx(kscattop+1,1)+D_H(kscattop+1)*Eyinc(kscattop+1,1) ; 
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Hy (kscattop+1,1)=Hy(kscattop+1, 1)- 
D_H(kscattop+1)*Exinc(kscattop+1,1); 

Hxinc(:,1)=C_H.*Hxinc(:,1)+D_H.*(Eyinc(zh, 1)-Eyinc(zh-1,1)) ; 

Hyinc(:,1)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)); 

end 

Dx (ze, 1) =C_E (ze) . *Dx (ze, 1) -D_E (ze) . * (Hy (ze, 1) -Hy (ze-1,1) ) ; 

Dy(ze,1)=C_E(ze) .*Dy(ze,1)+D_E(ze) .*(Hx(ze, 1)-Hx(ze-1,1)) ; 

Dx(kscatbot + 1,1)=Dx(kscatbot + 1,1)+D_E(kscatbot + 1)*Hyinc(kscatbot, 1) ; 

Dy(kscatbot+1,1)=Dy(kscatbot+1,1)-D_E(kscatbot+1)*Hxinc(kscatbot,1); 

Dx(kscattop+1,1)=Dx(kscattop+1,1)-D_E(kscattop+1)*Hyinc(kscattop+1,1) ; 
Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1, 1) ; 
Ex(ze,1)=epsinvxx(ze).*Dx(ze,1)+epsinvxy(ze).*Dy(ze,1)+epsinvxz(ze).*Dz 
(ze,1) ; 

Ey(ze,1)=epsinvxy(ze).*Dx(ze,1)+epsinvyy(ze).*Dy(ze,1)+epsinvyz(ze).*Dz 
(ze,1); 

Ez(ze, 1)=epsinvxz(ze) .*Dx(ze, 1)+epsinvyz (ze) .*Dy(ze, 1)+epsinvzz (ze) . *Dz 
(ze,1); 

Exinc(ze,1)=C_E (ze) .*Exinc(ze,1)-D_E(ze) .*(Hyinc(ze,1)-Hyinc (ze-1,1)); 
Eyinc(ze, 1)=C_E(ze) .*Eyinc(ze,1)+D_E(ze) . *(Hxinc(ze,1)-Hxinc(ze-1,1)); 
Exinc(kincstart+1,1)=-E0*sin(omega*Tbreak*dt)*(1-exp(-(Tbreak*dt- 
tdelay)^2/sigmadelay)); 

Hx(:,2)=C_H.*Hx(:,1)+D_H.*(Ey(zh,1)-Ey(zh-1,1)); 

Hy(:,2)=C_H.*Hy(:,1)-D_H.*(Ex(zh, 1)-Ex(zh-1,1)) ; 

Hx(kscatbot,2)=Hx(kscatbot,2)-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,2)=Hy(kscatbot,2)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,2)=Hx(kscattop+1,2)+D_H(kscattop+1)*Eyinc(kscattop+1,1) ; 
Hy(kscattop+1,2)=Hy(kscattop+1,2)-D_H(kscattop+1)*Exinc(kscattop+1,1) ; 
Hxinc(:,2)=C_H.*Hxinc(:,1)+D_H. * (Eyinc(zh,1)-Eyinc(zh-1,1)); 

Hyinc(:,2)=C_H.*Hyinc(:,1)-D_H. * (Exinc(zh,1)-Exinc(zh-1,1)); 

Hxinctilde(: , 1) = (Hxinc(ze, 1)+Hxinc(ze, 2)+Hxinc(ze-1,1)+Hxinc(ze- 

1.2) ) /4 ; 

Hyinctilde(:,1) = (Hyinc(ze,1)+Hyinc(ze, 2) +Hyinc(ze-1, 1)+Hyinc(ze- 

1.2) )/4; 

Hxtilde(:,1) = (Hx(ze, 1)+Hx(ze, 2) +Hx(ze-1, 1)+Hx(ze-1,2))/4 ; 

Hytilde(:,1) = (Hy(ze,1)+Hy(ze,2)+Hy(ze-1,1)+Hy (ze-l,2))/4; 

Poyntingx(: , 1)=-Ez(ze, 1) .*Hytilde(:,1); 

Poyntingy(:,l)=Ez(ze,l).*Hxtilde(:,1); 

Poyntingz(:,1)=Ex(ze,1).*Hytilde(:,1)-Ey(ze,1).*Hxtilde(:,1); 
PoyntingMag(:,1)=sqrt(Poyntingx(: , 1) . A 2+Poyntingy(: , 1) .^2+Poyntingz(:,1 
). A 2); 

for n=l:2*timeperiod 

Dx(ze, n+1)=C_E(ze) .*Dx(ze, n)-D_E(ze) .*(Hy(ze,n+1)-Hy(ze-l,n+l)); 

Dy(ze,n+1)=C_E(ze).*Dy(ze,n)+D_E(ze). *(Hx(ze,n+1)-Hx(ze-1,n+1)); 

Dx(kscatbot+1,n+1)=Dx(kscatbot+1,n+1)+D_E(kscatbot+1)*Hyinc(kscatbot,n+ 

i); 

Dy(kscatbot + 1,n+1)=Dy(kscatbot + 1, n+1)- 
D_E(kscatbot + 1)*Hxinc(kscatbot, n+1) ; 

Dx(kscattop+1,n+1)=Dx(kscattop+1, n+1)- 
D_E(kscattop+1)*Hyinc(kscattop+1, n+1) ; 

Dy(kscattop+1,n+1)=Dy(kscattop+1,n+1)+D_E(kscattop+1)*Hxinc(kscattop+1, 
n+1) ; 
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Ex(ze,n+1)=epsinvxx(ze).*Dx(ze,n+1)+epsinvxy(ze).*Dy(ze,n+1)+epsinvxz(z 
e).*Dz(ze,n+1); 

Ey(ze, n+1)=epsinvxy(ze) .*Dx(ze, n+1)+epsinvyy(ze) .*Dy(ze, n+1)+epsinvyz (z 
e).*Dz(ze,n+1); 

Ez(ze, n+1)=epsinvxz(ze) .*Dx(ze, n+1)+epsinvyz(ze) .*Dy(ze,n+1)+epsinvzz (z 
e).*Dz(ze,n+1); 

Exinc(ze,n+1)=C_E(ze) .*Exinc(ze,n)-D_E(ze) . *(Hyinc(ze,n+1)-Hyinc (ze- 
1,n+1)); 

Eyinc(ze,n+1)=C_E(ze).*Eyinc(ze,n)+D_E(ze). *(Hxinc(ze,n+1)-Hxinc(ze- 
1,n+1)); 

Exinc(kincstart+1,n+1)=-EO*sin(omega*(Tbreak+n)*dt)*(1-exp(- 
( (Tbreak+n) *dt-tdelay) A 2/sigmadelay) ) ; 

Hx(:,n+2)=C_H.*Hx(:,n+1)+D_H.*(Ey(zh,n+1)-Ey(zh-1, n+1)); 

Hy(:,n+2)=C_H.*Hy(:,n+1)-D_H.*(Ex(zh, n+1)-Ex(zh-1, n+1)) ; 

Hx(kscatbot,n+2)=Hx(kscatbot, n+2)- 
D_H(kscatbot)*Eyinc(kscatbot+ 1, n+1) ; 

Hy(kscatbot,n+2)=Hy(kscatbot,n+2)+D_H(kscatbot)*Exinc(kscatbot+1,n+1); 

Hx(kscattop+1,n+2)=Hx(kscattop+1,n+2)+D_H(kscattop+1)*Eyinc(kscattop+1, 
n+1) ; 

Hy (kscattop+1,n+2)=Hy(kscattop+1, n+2)- 
D_H(kscattop+1)*Exinc(kscattop+1, n+1); 

Hxinc(:,n+2)=C_H.*Hxinc(:,n+1)+D_H.*(Eyinc(zh, n+1)-Eyinc(zh-1, n+1)) ; 
Hyinc(:,n+2)=C_H.*Hyinc(:,n+1)-D_H.*(Exinc(zh,n+1)-Exinc(zh-1,n+1)); 
Hxinctilde(:,n+1)=(Hxinc(ze,n+1)+Hxinc(ze,n+2)+Hxinc(ze- 
1,n+1)+Hxinc(ze-l,n+2))/4; 

Hyinctilde(: , n+1) = (Hyinc(ze,n+1)+Hyinc(ze,n+2)+Hyinc(ze- 
1,n+1)+Hyinc(ze-l,n+2))/4; 

Hxtilde(: , n+1) = (Hx(ze, n+1)+Hx(ze, n+2)+Hx(ze-1, n+1)+Hx(ze-1, n+2))/4 ; 
Hytilde(:,n+1)=(Hy(ze,n+1)+Hy(ze,n+2)+Hy(ze-1,n+1)+Hy(ze-1,n+2))/4; 
Poyntingx(:,n+1)=-Ez(ze,n+1).*Hytilde(:,n+1); 

Poyntingy(:,n+1)=Ez(ze,n+1).*Hxtilde(:,n+1); 

Poyntingz(:,n+1)=Ex(ze,n+1).*Hytilde(:,n+1)- 
Ey(ze,n+1).*Hxtilde(:,n+1); 

PoyntingMag(: , n+1)=sqrt(Poyntingx(: , n+1) . A 2+Poyntingy(:,n+1) . A 2+Poyntin 

gz(:,n+1) . A 2); 

end 

Intensity=PoyntingMag(: , 1) ; 
for n=2:2*timeperiod 

Intensity=Intensity+2*PoyntingMag(:,n); 

end 

Intensity=Intensity+PoyntingMag(:,2*timeperiod+l); 
Intensity=Intensity*dt/EO A 2/(4*pi/omega); 

kPMLbot3=kPMLbot; kscatbot3=kscatbot; kairbot3=kairbot; 
kglassbot3=kglassbot; 

kglasstop3=kglasstop; kairtop3=kairtop; kscattop3=kscattop; 
kPMLtop3=kPMLtop; 

Ntot3=Ntot; dz3=dz; dt3=dt; 

Ex3=Ex; Ey3=Ey; Ez3=Ez; 
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Hx3=Hx; Hy3=Hy; 

Exinc3=Exinc; Eyinc3=Eyinc; 

Hxinc3=Hxinc; Hyinc3=Hyinc; 

Hxinctilde3=Hxinctilde; Hyinctilde3=Hyinctilde; 

Hxtilde3=Hxtilde; Hytilde3=Hytilde; 

Poyntingx3=Poyntingx; Poyntingy3=Poyntingy; Poyntingz3=Poyntingz; 

PoynitngMag3=PoyntingMag; Intensity3=Intensity; 

timeperiod3=timeperiod; Tfin3=Tfin; Tbreak3=Tbreak; 

kincstart3=kincstart; 

zTotE3=dz*((1rNtottl)-kglassbot-1); 

zTotH3=dz*((1:Ntot)-1/2-kglassbot); 

toe 


tic 

%% Finer grid 
refinement=2; 

NLCP=NLCP*refinement %number of numerical cells in the LCP 

dz=h/NLCP; 


Nglass=Nglass*refinement; 

%number 

of 

numerical 

cells 

in 

the 

glass 

Nair=Nair*refinement; 

%number 

of 

numerical 

cells 

in 

the 

air 

between glass and scattering 
Nscat=Nscat*refinement; 

layer 

%number 

of 

numerical 

cells 

in 

the 


scattering later 
NPML=NPML*refinement; 

%number 

of 

numerical 

cells 

in 

the 



perfectly matched layer 
cdtoverdz = l; 
dt=dz*cdtoverdz; 

Tfin=floor(EndTime/dt)+1; 
timeperiod=round(2*pi/omega/dt); 

sigmamax=2/dt; %10 A 17; %in the PML, sigmaz=sigmamax((kPML- 

1)/(NPML-1)) A sigmam 
sigmam=3; 

kPMLbot=NPML; % sigma defined but equal to zero at 

k=kPMLbot 


kscatbot=kPMLbottNscat; 
total field; kscatbot-1/2 
kairbot=kscatbottNair; 
bottom plate 

kglassbot=kairbottNglass; 
plate/where the anchoring 
kglasstop=kglassbottNLCP; 
plate/where the anchoring 
kairtop=kglasstoptNglass; 
plate 

kscattop=kairtoptNair; 
kscattoptl/2 is not 
kPMLtop=kscattoptNscat; 
k=kPMLtop 
Ntot=kPMLtoptNPML; 
ze=2:Ntot; 

in the D and E equations 
zh=2:Ntottl; 

z4 = linspace(0,h,NLCPtl) ; 
kincstart=refinement*(kinc 


% k=kscatbot is considered to be in the 
is in the scattered field 

% k=kairbot is the bottom edge of the 

% k=kglassbot is the top edge of the bottom 
condition is applied 

% k=kglasstop is the bottom edge of the top 
condition is applied 

% k=kairtop is the top edge of the top 

% k=kscattop is in the total field; 

% sigma defined but equal to zero at 

% Ntot cells in total% ze=2:Ntot; 

% k values for updating the gradients of H 


%z values used for computing the orientation 
start3-kglassbot3)tkglassbot; 
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%% compute the permittivity 
solinit=bvpinit(z4, @LEodeinit); 
sol=bvp4c(SLEodefun,@LEbcfun, solinit) ; 
u=deval(sol,z4) ; 
psiLE=u (1, : ); 
thetaLE=u (3, : ); 

nl=cos(thetaLE).*cos(psiLE); 
n2=cos(thetaLE).*sin(psiLE); 
n3=sin(thetaLE); 

epsxx=zeros(Ntot + 1, 1) ; 
epsyy=zeros(Ntot + 1,1) ; 
epszz=zeros(Ntot+1,1); 
epsxy=zeros(Ntot + 1, 1) ; 
epsxz=zeros(Ntot + 1,1) ; 
epsyz=zeros(Ntot+1,1); 

epsxx(1:kPMLbot+1)=1; epsyy(1:kPMLbot+1)=1; epszz(1:kPMLbot+1)=1; 

epsxx(kPMLbot+2:kairbot)=epsilon_air; 

epsyy(kPMLbot+2:kairbot)=epsilon_air; 

epszz(kPMLbot+2:kairbot)=epsilon_air; 

epsxx(kairbot + 1:kglassbot)=epsilon_glass; 

epsyy(kairbot + 1:kglasshot)=epsilon_glass; 

epszz (kairbot + 1:kglassbot)=epsilon_glass; 

epsxx(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*nl. A 2; 

epsyy(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n2. A 2; 

epszz(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n3. A 2; 

epsxy(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n2; 

epsxz(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n3; 

epsyz(kglassbot+1:kglasstop+1)=delta_epsilon*n2.*n3; 

epsxx(kglasstop+2:kairtop+1)=epsilon_glass; 

epsyy(kglasstop+2:kairtop+1)=epsilon_glass; 

epszz(kglasstop+2:kairtop+1)=epsilon_glass; 

epsxx(kairtop+2:kPMLtop)=epsilon_air; 

epsyy(kairtop+2:kPMLtop)=epsilon_air; 

epszz(kairtop+2:kPMLtop)=epsilon_air; 

epsxx(kPMLtop+1:Ntot)=1; epsyy(kPMLtop+1:Ntot)=1; 

epszz(kPMLtop+1:Ntot)=1; 

epsdet=epsxx.*(epsyy.*epszz-epsyz. A 2)-epsxy.*(epsxy.*epszz- 

epsxz.*epsyz)+epsxz.*(epsxy.*epsyz-epsyy.*epsxz); 

epsinvxx=(epsyy.*epszz-epsyz. A 2)./epsdet; 

epsinvxy=(epsxz.*epsyz-epsxy.*epszz)./epsdet; 

epsinvxz=(epsxy.*epsyz-epsyy.*epsxz)./epsdet; 

epsinvyy=(epsxx.*epszz-epsxz. A 2)./epsdet; 

epsinvyz=(epsxy.*epsxz-epsxx.*epsyz)./epsdet; 

epsinvzz=(epsxx.*epsyy-epsxy. A 2)./epsdet; 

%% Define coefficients 

C_E=ones(Ntot + 1,1); % factor of the D A n term in the D 

update equations 

D_E=cdtoverdz*ones(Ntot+1,1); % factor of the curl H term in the D 

update equations 
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factor of the H A n term in the H 


C_H=ones(Ntot,1); 
update equations 

D_H=cdtoverdz*ones(Ntot, 1) ; % factor of the curl H term in the H 

update equations 

for k=l:kPMLbot % fill in the effect of sigma in the 

PML 

sigmaz=sigmamax*((kPMLbot-k+1)/kPMLbot) A sigmam; 

C_E(k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(k)=2*cdtoverdz/(2+dt*sigmaz) ; 

C_E(Ntot+2-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(Ntot+2-k)=2*cdtoverdz/(2+dt*sigmaz); 
sigmaz=sigmamax*((kPMLbot-k+0.5)/kPMLbot) A sigmam; 

C_H(k) = (2-dt*sigmaz)/(2+dt*sigmaz) ; 

D_H(k)=2*cdtoverdz/(2+dt*sigmaz); 

C_H(Ntot+l-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_H(Ntot+l-k)=2*cdtoverdz/(2+dt*sigmaz); 

end 


Hz is zero and not stored. 


%% Define electromagnetic field variables 
Dx=zeros(Ntot + 1,2*timeperiod+l) ; % 

time grid 

Dy=zeros(Ntot + 1,2*timeperiod+l) ; % 

initially (n=0), but this is not stored. 

Dz=zeros(Ntot + 1,2*timeperiod+l) ; % 

k=Ntot+l but they are stored so that 
Ex=zeros(Ntot + 1,2*timeperiod+l) ; 
z=(k-kglassbot-1)*dz, 

Ey=zeros(Ntot + 1,2*timeperiod+l) ; 

Ex(k+1,n) 

Ez=zeros(Ntot + 1,2*timeperiod+l) ; 

Hx=zeros(Ntot,2*timeperiod+2 ) ; 
z=(k-kglassbot+1/2)*dz 
Hy=zeros(Ntot,2*timeperiod+2 ) ; 

Hx(k+1,n+1) 

Exinc=zeros(Ntot+1,2*timeperiod+l); 
can be nonzero, depending on epsilon 
Eyinc=zeros(Ntot+1,2*timeperiod+l); 

Hxinc=zeros(Ntot,2*timeperiod+2 ) ; 
are the incident field 
Hyinc=zeros(Ntot,2*timeperiod+2) ; 

Hxinctilde=zeros(Ntot-1,2*timeperiod+l); 

Hxinc onto the grid 

Hyinctilde=zeros(Ntot-1,2*timeperiod+l) ; 

Hxtilde=zeros(Ntot-1,2*timeperiod+l); 
onto the grid 

Hytilde=zeros(Ntot-1,2*timeperiod+l); 

Poyntingx=zeros(Ntot-1,2*timeperiod+l) ; 
on the space-time grid, and so H must be 
Poyntingy=zeros(Ntot-1,2*timeperiod+l); 
space. Thus it is not defined at k=l or k=Ntot+l 
Poyntingz=zeros(Ntot-1,2*timeperiod+l); % or at t=Tfin 

PoyntingMag=zeros(Ntot-1,2*timeperiod+l); 

Tbreak=Tfin-2*timeperiod-l; 


D and E are known on the space- 

D and E fields are zero 

Dx and Ex are zero at k=l and 

Dx(k,n) is Dx at t=n*dt and 

or Ex_k A n is stored as 

H is known on the half-grid. 
Hx(k,n) is Hx at t=(n-1/2)*dt, 

Hx_(k+1/2) A (n+1/2 ) is stored as 


Ez 


Exinc, Eyinc, Hxinc, and Hyinc 

% Hxinctilde interpolates 

% Hxtilde interpolates Hx 

% The Poynting vector is known 
interpolated in time and in 


The starting main loop, not storing every step 
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Exinc(kincstart + 1,1)=-EO*sin(omega*dt)*(1-exp(-(dt- 
tdelay) A 2/sigmadelay)); 

Hx(kscatbot,1)=-D_H(kscatbot)*Eyinc(kscatbot+1,1); 

Hy(kscatbot,1)=D_H(kscatbot)*Exinc(kscatbot + 1,1); 

Hxinc(:,1)=D_H .* (Eyinc(zh, 1)-Eyinc(zh-1,1) ) ; 

Hyinc(:,1)=-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)) ; 
for n=2:Tbreak-l 

Dx(ze,1)=C_E(ze) .*Dx(ze,1)-D_E(ze) .*(Hy(ze, 1)-Hy(ze-1,1)) ; 

Dy (ze, 1) =C_E (ze) . *Dy (ze, 1) +D_E (ze) . * (Hx (ze, 1) -Hx (ze-1, 1) ) ; 

Dx(kscatbot+1,1)=Dx(kscatbot+1,1)+D_E(kscatbot+1)*Hyinc(kscatbot,1); 
Dy(kscatbot+1,1)=Dy(kscatbot+1,1)-D_E(kscatbot+1)*Hxinc(kscatbot,1); 
Dx(kscattop+1,1)=Dx(kscattop+1,1)- 
D_E(kscattop+1)*Hyinc(kscattop+1, 1) ; 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1); 

Ex(ze, 1)=epsinvxx(ze) .*Dx(ze, 1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz(ze) .*Dz 
(ze,1) ; 

Ey(ze,1)=epsinvxy(ze).*Dx(ze,1)+epsinvyy(ze).*Dy(ze,1)+epsinvyz(ze).*Dz 
(ze,1) ; 

Exinc(ze,1)=C_E(ze).*Exinc(ze,1)-D_E(ze).*(Hyinc(ze,1)-Hyinc(ze- 

i,D); 

Eyinc(ze,1)=C_E(ze) .*Eyinc(ze, 1)+D_E(ze) .*(Hxinc(ze,1)-Hxinc (ze- 

i,D); 

Exinc(kincstart+1,1)=-E0*sin(omega*n*dt)*(1-exp(-(n*dt- 
tdelay) A 2/sigmadelay)) ; 

Hx(:,1)=C_H.*Hx(:,1)+D_H.*(Ey(zh,1)-Ey(zh-1,1)); 

Hy(:,1)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1, 1)) ; 

Hx(kscatbot,1)=Hx(kscatbot,1)-D_H(kscatbot)*Eyinc(kscatbot + 1,1) ; 

Hy(kscatbot,1)=Hy(kscatbot,1)+D_H(kscatbot)*Exinc(kscatbot + 1, 1) ; 

Hx(kscattop+1,1)=Hx(kscattop+1,1)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 

Hy(kscattop+1,1)=Hy(kscattop+1, 1)- 
D_H(kscattop+1)*Exinc(kscattop+1,1); 

Hxinc(:,1)=C_H.*Hxinc(:,1)+D_H.*(Eyinc(zh, 1)-Eyinc(zh-1,1)) ; 

Hyinc(:,1)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh, 1)-Exinc(zh-1,1)) ; 

end 

%n=Tbreak; 

Dx(ze,1)=C_E(ze) .*Dx(ze,1)-D_E(ze) .* (Hy (ze,1)-Hy(ze-1,1)); 

Dy (ze, 1) =C_E (ze) . *Dy (ze, 1) +D_E (ze) . * (Hx (ze, 1) -Hx (ze-1,1) ) ; 

Dx(kscatbot + 1,1)=Dx(kscatbot + 1,1)+D_E(kscatbot + 1)*Hyinc(kscatbot, 1) ; 

Dy(kscatbot+1,1)=Dy(kscatbot+1,1)-D_E(kscatbot+1)*Hxinc(kscatbot,1); 

Dx(kscattop+1,1)=Dx(kscattop+1,1)-D_E(kscattop+1)*Hyinc(kscattop+1,1); 
Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1) ; 
Ex(ze,1)=epsinvxx(ze) .*Dx(ze,1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz (ze) .*Dz 
(ze,1); 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze, 1)+epsinvyy (ze) .*Dy(ze, 1)+epsinvyz (ze) . *Dz 
(ze,1) ; 

Ez(ze,1)=epsinvxz(ze).*Dx(ze,1)+epsinvyz(ze).*Dy(ze,1)+epsinvzz(ze).*Dz 
(ze,1) ; 

Exinc(ze,1)=C_E(ze) .*Exinc(ze,1)-D_E(ze) .*(Hyinc(ze, 1)-Hyinc(ze-1, 1)) ; 
Eyinc(ze, 1)=C_E(ze) .*Eyinc(ze,1)+D_E (ze) .*(Hxinc(ze,1)-Hxinc(ze-1,1)); 
Exinc(kincstart + 1,1)=-E0*sin(omega*Tbreak*dt)*(1-exp(-(Tbreak*dt- 
tdelay) A 2/sigmadelay)); 
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Hx ( : , 2) =C_H. *Hx ( : , 1) +D_H. * (Ey (zh, 1) -Ey (zh-1, 1) ) ; 

Hy(:,2)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1,1)); 

Hx(kscatbot,2)=Hx(kscatbot,2)-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,2)=Hy(kscatbot,2)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,2)=Hx(kscattop+1,2)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 
Hy(kscattop+1,2)=Hy(kscattop+1,2)-D_H(kscattop+1)*Exinc(kscattop+1, 1) ; 
Hxinc(:,2)=C_H.*Hxinc(:,1)+D_H.* (Eyinc(zh, 1)-Eyinc(zh-1, 1)) ; 

Hyinc(:,2)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)); 

Hxinctilde(:,1)=(Hxinc(ze,1)+Hxinc(ze,2)+Hxinc(ze-1,1)+Hxinc(ze- 
1,2))/4; 

Hyinctilde(: , 1) = (Hyinc(ze,1)+Hyinc(ze,2)+Hyinc(ze-1,1)+Hyinc (ze- 
1,2))/4; 

Hxtilde(: , 1) = (Hx(ze,1)+Hx(ze, 2 ) +Hx(ze-1,1)+Hx(ze-l,2))/4; 

Hytilde(:,l) = (Hy(ze,l)+Hy(ze,2)+Hy(ze-1,1)+Hy (ze-1,2))/4; 

Poyntingx(: , 1)=-Ez(ze, 1) .*Hytilde(:,1); 

Poyntingy(:,1)=Ez(ze, 1) .*Hxtilde(:,1); 

Poyntingz(: , 1)=Ex(ze, 1) .*Hytilde(: , 1)-Ey(ze, 1) .*Hxtilde(:, 1) ; 
PoyntingMag(: , 1)=sqrt(Poyntingx(: , 1) . A 2+Poyntingy(:,1) . A 2+Poyntingz(:,1 
). A 2); 

%for n=Tbreak:Tfin 

%makes Dx(:,l) actually at n=Tfin-l-2*timeperiod 
for n=l:2*timeperiod 

Dx(ze,n+1)=C_E(ze).*Dx(ze,n)-D_E(ze).*(Hy(ze,n+1)-Hy(ze-l,n+l)); 

Dy(ze,n+1)=C_E(ze).*Dy(ze,n)+D_E(ze). *(Hx(ze,n+1)-Hx(ze-1,n+1)); 

Dx(kscatbot + 1, n+1)=Dx(kscatbot + 1,n+1)+D_E(kscatbot + 1)*Hyinc(kscatbot,n+ 

i); 

Dy (kscatbot + 1,n+1)=Dy(kscatbot + 1, n+1)- 
D_E(kscatbot+1)*Hxinc(kscatbot,n+1); 

Dx(kscattop+1,n+1)=Dx(kscattop+1, n+1)- 
D_E(kscattop+1)*Hyinc(kscattop+1, n+1) ; 

Dy(kscattop+1, n+1)=Dy(kscattop+1,n+1)+D_E(kscattop+1)*Hxinc(kscattop+1, 
n+1) ; 

Ex(ze,n+1)=epsinvxx(ze) .*Dx(ze, n+1)+epsinvxy(ze) .*Dy(ze, n+1)+epsinvxz(z 
e).*Dz(ze,n+1); 

Ey(ze,n+1)=epsinvxy(ze).*Dx(ze,n+1)+epsinvyy(ze).*Dy(ze,n+1)+epsinvyz(z 
e) . *Dz (ze,n+1); 

Ez(ze, n+1)=epsinvxz(ze) .*Dx(ze, n+1)+epsinvyz(ze) .*Dy(ze,n+1)+epsinvzz(z 
e).*Dz(ze,n+1); 

Exinc(ze,n+1)=C_E (ze) .*Exinc(ze,n)-D_E(ze) . *(Hyinc(ze,n+1)-Hyinc (ze- 
1,n+1)); 

Eyinc(ze, n+1)=C_E(ze) .*Eyinc(ze,n)+D_E(ze) .*(Hxinc(ze,n+1)-Hxinc (ze- 
1,n+1)); 

Exinc(kincstart + 1, n+1)=-E0*sin(omega*(Tbreak+n)*dt)*(1-exp(- 
( (Tbreak+n)*dt-tdelay) A 2/sigmadelay)); 

Hx(:,n+2)=C_H.*Hx(:,n+1)+D_H.*(Ey(zh, n+1)-Ey(zh-1, n+1)) ; 

Hy(:,n+2)=C_H.*Hy(:,n+1)-D_H.*(Ex(zh,n+1)-Ex(zh-1, n+1)) ; 

Hx(kscatbot,n+2)=Hx(kscatbot, n+2)- 
D_H(kscatbot)*Eyinc(kscatbot+1,n+1); 

Hy(kscatbot,n+2)=Hy(kscatbot,n+2)+D_H(kscatbot)*Exinc(kscatbot + 1,n+1) ; 
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Hx(kscattop+1, n+2)=Hx(kscattop+1,n+2)+D_H(kscattop+1)*Eyinc(kscattop+1, 
n+1) ; 

Hy(kscattop+1,n+2)=Hy(kscattop+1, n+2)- 
D_H(kscattop+1)*Exinc(kscattop+1, n+1) ; 

Hxinc(:,n+2)=C_H.*Hxinc(:,n+1)+D_H.*(Eyinc(zh,n+1)-Eyinc(zh-1, n+1)) ; 
Hyinc(:,n+2)=C_H.*Hyinc(:,n+1)-D_H.*(Exinc(zh,n+1)-Exinc(zh-1, n+1)) ; 
Hxinctilde(: , n+1) = (Hxinc(ze,n+1)+Hxinc(ze,n+2)+Hxinc (ze- 
1,n+1)+Hxinc(ze-l,n+2))/4; 

Hyinctilde(:,n+1) = (Hyinc(ze,n+1)+Hyinc(ze,n+2)+Hyinc (ze- 
1,n+1)+Hyinc(ze-l,n+2))/4; 

Hxtilde(:,n+1) = (Hx(ze,n+1)+Hx(ze,n+2)+Hx(ze-1,n+1)+Hx(ze-1, n+2))/4 ; 
Hytilde(: , n+1) = (Hy(ze,n+1)+Hy(ze,n+2)+Hy(ze-1,n+1)+Hy(ze-1,n+2))/4; 
Poyntingx(:,n+1)=-Ez(ze,n+1).*Hytilde(:,n+1); 

Poyntingy(:,n+1)=Ez(ze,n+1).*Hxtilde(:,n+1); 

Poyntingz(:,n+1)=Ex(ze,n+1).*Hytilde(:,n+1)- 
Ey(ze,n+1).*Hxtilde(:,n+1); 

PoyntingMag(:,n+1)=sqrt(Poyntingx(:,n+1). A 2+Poyntingy(:,n+1). A 2+Poyntin 

gz(:,n+1) . A 2); 

end 

Intensity=PoyntingMag(: , 1) ; 
for n=2:2*timeperiod 

Intensity=Intensity+2*PoyntingMag(: , n) ; 

end 

Intensity=Intensity+PoyntingMag(:,2*timeperiod+l); 
Intensity=Intensity*dt/EO A 2/(4*pi/omega); 

kPMLbot4=kPMLbot; kscatbot4=kscatbot; kairbot4=kairbot; 
kglassbot4=kglassbot; 

kglasstop4=kglasstop; kairtop4=kairtop; kscattop4=kscattop; 
kPMLtop4=kPMLtop; 

Ntot4=Ntot; dz4=dz; dt4=dt; 

Ex4=Ex; Ey4=Ey; Ez4=Ez; 

Hx4=Hx; Hy4=Hy; 

Exinc4=Exinc; Eyinc4=Eyinc; 

Hxinc4=Hxinc; Hyinc4=Hyinc; 

Hxinctilde4=Hxinctilde; Hyinctilde4=Hyinctilde; 

Hxtilde4=Hxtilde; Hytilde4=Hytilde; 

Poyntingx4=Poyntingx; Poyntingy4=Poyntingy; Poyntingz4=Poyntingz; 

PoynitngMag4=PoyntingMag; Intensity4=Intensity; 

timeperiod4=timeperiod; Tfin4=Tfin; Tbreak4=Tbreak; 

kincstart4=kincstart; 

zTotE4=dz*((1:Ntot+l)-kglassbot-1); 

zTotH4=dz*((1:Ntot)-1/2-kglassbot); 

toe 

tic 

%% Finer grid 
refinement=2; 

NLCP=NLCP*refinement %number of numerical cells in the LCP 

dz=h/NLCP; 

Nglass=Nglass*refinement; %number of numerical cells in the glass 

Nair=Nair*refinement; %number of numerical cells in the air 

between glass and scattering layer 
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%number of numerical cells in the 


Nscat=Nscat*refinement; 
scattering later 

NPML=NPML*refinement; %number of numerical cells in the 

perfectly matched layer 
cdtoverdz = l; 
dt=dz*cdtoverdz; 

Tfin=floor(EndTime/dt)+1; 
timeperiod=round(2*pi/omega/dt) ; 

sigmamax=2/dt; %10 A 17; %in the PML, sigmaz = sigmamax ( (kPML- 

1)/(NPML-1)) A sigmam 
sigmam=3; 

kPMLbot=NPML; % sigma defined but equal to zero at 

k=kPMLbot 


kscatbot=kPMLbot+Nscat; 
total field; kscatbot-1/2 
kairbot=kscatbot+Nair; 
bottom plate 

kglassbot=kairbot+Nglass; 
plate/where the anchoring 
kglasstop=kglassbot+NLCP; 
plate/where the anchoring 
kairtop=kglasstop+Nglass; 
plate 

kscattop=kairtop+Nair; 
kscattop+1/2 is not 
kPMLtop=kscattop+Nscat; 
k=kPMLtop 
Ntot=kPMLtop+NPML; 
ze=2:Ntot; 

in the D and E equations 
zh=2:Ntot+1; 

z5=linspace(0,h,NLCP+1); 
kincstart=refinement*(kinc 


% k=kscatbot is considered to be in the 
is in the scattered field 

% k=kairbot is the bottom edge of the 

% k=kglassbot is the top edge of the bottom 
condition is applied 

% k=kglasstop is the bottom edge of the top 
condition is applied 

% k=kairtop is the top edge of the top 

% k=kscattop is in the total field; 

% sigma defined but equal to zero at 

% Ntot cells in total% ze=2:Ntot; 

% k values for updating the gradients of H 


%z values used for computing the orientation 
start4-kglassbot4)+kglassbot; 


%% compute the permittivity 
solinit=bvpinit(z5,@LEodeinit); 
sol=bvp4c(@LEodefun,@LEbcfun,solinit); 
u=deval(sol,z5) ; 
psiLE=u (1, : ); 
thetaLE=u (3, : ); 


nl=cos(thetaLE).*cos(psiLE); 
n2=cos(thetaLE).*sin(psiLE); 
n3=sin(thetaLE); 


epsxx=zeros(Ntot + 1, 1) ; 
epsyy=zeros(Ntot + 1,1) ; 
epszz=zeros(Ntot+1,1); 
epsxy=zeros(Ntot + 1, 1) ; 
epsxz=zeros(Ntot+1,1); 
epsyz=zeros(Ntot + 1, 1) ; 

epsxx(1:kPMLbot+1)=1; epsyy(1:kPMLbot+1)=1; epszz(1:kPMLbot+1)=1; 
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epsxx(kPMLbot+2:kairbot)=epsilon_air; 
epsyy(kPMLbot+2:kairbot)=epsilon_air; 
epszz(kPMLbot+2:kairbot)=epsilon_air; 
epsxx(kairbot+1:kglassbot)=epsilon_glass; 
epsyy(kairbot+ 1:kglassbot)=epsilon_glass; 
epszz(kairbot + 1:kglassbot)=epsilon_glass; 

epsxx(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*nl. A 2; 

epsyy(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n2. A 2; 

epszz(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n3. A 2; 

epsxy(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n2; 

epsxz(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n3; 

epsyz(kglassbot+1:kglasstop+1)=delta_epsilon*n2.*n3; 

epsxx(kglasstop+2:kairtop+1)=epsilon_glass; 

epsyy(kglasstop+2:kairtop+1)=epsilon_glass; 

epszz(kglasstop+2:kairtop+1)=epsilon_glass; 

epsxx(kairtop+2:kPMLtop)=epsilon_air; 

epsyy(kairtop+2:kPMLtop)=epsilon_air; 

epszz(kairtop+2:kPMLtop)=epsilon_air; 

epsxx(kPMLtop+1:Ntot)=1; epsyy(kPMLtop+1:Ntot)=1; 

epszz(kPMLtop+1:Ntot)=1; 


epsdet=epsxx.*(epsyy.*epszz-epsyz. A 2)-epsxy.*(epsxy.*epszz- 

epsxz.*epsyz)+epsxz.*(epsxy.*epsyz-epsyy.*epsxz); 

epsinvxx=(epsyy.*epszz-epsyz. A 2)./epsdet; 

epsinvxy=(epsxz.*epsyz-epsxy.*epszz)./epsdet; 

epsinvxz=(epsxy.*epsyz-epsyy.*epsxz)./epsdet; 

epsinvyy=(epsxx.*epszz-epsxz. A 2)./epsdet; 

epsinvyz=(epsxy.*epsxz-epsxx.*epsyz)./epsdet; 

epsinvzz=(epsxx.*epsyy-epsxy. A 2)./epsdet; 


%% Define coefficients 

C_E=ones(Ntot + 1, 1); 

update equations 

D_E=cdtoverdz*ones(Ntot + 1,1) ; 

update equations 

C_H=ones(Ntot,1); 

update equations 

D_H=cdtoverdz*ones(Ntot, 1) ; 

update equations 

for k=l:kPMLbot 

PML 


% factor of the D A n term in the D 

% factor of the curl H term in the D 

% factor of the H A n term in the H 

% factor of the curl H term in the H 

% fill in the effect of sigma in the 


sigmaz=sigmamax*((kPMLbot-k+1)/kPMLbot) A sigmam; 
C_E(k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(k)=2*cdtoverdz/(2+dt*sigmaz); 

C_E(Ntot+2-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(Ntot+2-k)=2*cdtoverdz/(2+dt*sigmaz); 
sigmaz=sigmamax*((kPMLbot-k+0.5)/kPMLbot) A sigmam; 
C_H(k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_H(k)=2*cdtoverdz/(2+dt*sigmaz); 

C_H(Ntot+l-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_H(Ntot+l-k)=2*cdtoverdz/(2+dt*sigmaz); 

end 


%% Define electromagnetic field variables 
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Dx=zeros(Ntot + 1,2*timeperiod+l) ; % 

time grid 

Dy=zeros(Ntot + 1,2*timeperiod+l) ; % 

initially (n=0), but this is not stored. 
Dz=zeros(Ntot + 1,2*timeperiod+l) ; % 

k=Ntot+l but they are stored so that 
Ex=zeros(Ntot + 1,2*timeperiod+l) ; % 

z=(k-kglassbot-1)*dz, 

Ey=zeros(Ntot + 1,2*timeperiod+l) ; % 

Ex(k+1,n) 

Ez=zeros(Ntot+1,2*timeperiod+l); % 

Hx=zeros(Ntot,2*timeperiod+2 ) ; % 

z=(k-kglassbot+1/2)*dz 

Hy=zeros(Ntot,2*timeperiod+2 ) ; % 

Hx(k+1,n+1) 

Exinc=zeros(Ntot + 1,2*timeperiod+l) ; % 

can be nonzero, dependending on epsilon 
Eyinc=zeros(Ntot + 1,2*timeperiod+l) ; % 

Hxinc=zeros(Ntot,2*timeperiod+2); % 

are the incident field 
Hyinc=zeros(Ntot,2*timeperiod+2) ; 
Hxinctilde=zeros(Ntot-1,2*timeperiod+l) ; 
Hxinc onto the grid 

Hyinctilde=zeros(Ntot-1,2*timeperiod+l) ; 
Hxtilde=zeros(Ntot-1,2*timeperiod+l); 
onto the grid 

Hytilde=zeros(Ntot-1,2*timeperiod+l); 
Poyntingx=zeros(Ntot-1,2*timeperiod+l) ; 
on the space-time grid, and so H must be 
Poyntingy=zeros(Ntot-1,2*timeperiod+l) ; 
space. Thus it is not defined at k=l or 
Poyntingz = zeros(Ntot-1,2*timeperiod+l) ; 
PoyntingMag=zeros(Ntot-1,2*timeperiod+l) 
Tbreak=Tfin-2*timeperiod-l; 


D and E are known on the space- 

D and E fields are zero 

Dx and Ex are zero at k=l and 

Dx(k,n) is Dx at t=n*dt and 

or Ex_k A n is stored as 

H is known on the half-grid. 
Hx(k,n) is Hx at t=(n-1/2)*dt, 

Hx_(k+1/2) A (n+1/2) is stored as 

Hz is zero and not stored. Ez 

Exinc, Eyinc, Hxinc, and Hyinc 

% Hxinctilde interpolates 

% Hxtilde interpolates Hx 

% The Poynting vector is known 

% interpolated in time and in 

k=Ntot+l 

% or at t=Tfin 


% The starting main loop, not storing every step 
Exinc(kincstart + 1,1)=-E0*sin(omega*dt)*(1-exp(-(dt- 
tdelay) A 2/sigmadelay)); 

Hx(kscatbot,1)=-D_H(kscatbot)*Eyinc(kscatbot + 1,1) ; 

Hy(kscatbot,1)=D_H(kscatbot)*Exinc(kscatbot + 1,1) ; 

Hxinc(:,1)=D_H .* (Eyinc(zh,1)-Eyinc(zh-1, 1)) ; 

Hyinc(:,1)=-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)) ; 
for n=2:Tbreak-l 

Dx(ze,1)=C_E(ze).*Dx(ze,1)-D_E(ze).*(Hy(ze,1)-Hy(ze-1,1)); 

Dy (ze, 1) =C_E (ze) . *Dy (ze, 1) +D_E (ze) . * (Hx (ze, 1) -Hx (ze-1, 1) ) ; 

Dx(kscatbot + 1,1)=Dx(kscatbot + 1,1)+D_E(kscatbot + 1)*Hyinc(kscatbot, 1) ; 
Dy(kscatbot+1,1)=Dy(kscatbot+1,1)-D_E(kscatbot+1)*Hxinc(kscatbot,1); 
Dx (kscattop+1,1)=Dx(kscattop+1, 1)- 
D_E(kscattop+1)*Hyinc(kscattop+1,1) ; 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1); 

Ex(ze,1)=epsinvxx(ze) .*Dx(ze,1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz (ze) .*Dz 
(ze,1); 


91 


Ey(ze,1)=epsinvxy(ze).*Dx(ze,1)+epsinvyy(ze).*Dy(ze,1)+epsinvyz(ze).*Dz 
(ze,1) ; 

Exinc(ze,1)=C_E(ze) .*Exinc(ze, 1)-D_E(ze) .*(Hyinc(ze,1)-Hyinc(ze- 

i,D); 

Eyinc(ze,1)=C_E(ze) .*Eyinc(ze, 1)+D_E(ze) .*(Hxinc(ze,1)-Hxinc(ze- 

i,D); 

Exinc(kincstart + 1,1)=-EO*sin(omega*n*dt)*(1-exp(-(n*dt- 
tdelay)^2/sigmadelay)) ; 

Hx(:,1)=C_H.*Hx(:,1)+D_H.*(Ey(zh,1)-Ey(zh-1,1)); 

Hy ( : , 1) =C_H. *Hy ( : , 1) -D_H. * (Ex (zh, 1) -Ex (zh-1, 1) ) ; 

Hx(kscatbot,1)=Hx(kscatbot,1)-D_H(kscatbot)*Eyinc(kscatbot+1,1); 

Hy(kscatbot,1)=Hy(kscatbot,1)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,1)=Hx(kscattop+1,1)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 

Hy (kscattop+1,1)=Hy(kscattop+1, 1)- 
D_H(kscattop+1)*Exinc(kscattop+1,1) ; 

Hxinc(:,1)=C_H.*Hxinc(:,1)+D_H.*(Eyinc(zh,1)-Eyinc(zh-1,1)) ; 

Hyinc(:,1)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh, 1)-Exinc(zh-1,1)) ; 

end 

%n=Tbreak; 

Dx(ze,1)=C_E(ze) .*Dx(ze,1)-D_E(ze) .*(Hy(ze, 1)-Hy(ze-1,1)) ; 

Dy(ze,1)=C_E(ze).*Dy(ze,1)+D_E(ze).*(Hx(ze,1)-Hx(ze-1,1)); 

Dx(kscatbot+1,1)=Dx(kscatbot+1,1)+D_E(kscatbot+1)*Hyinc(kscatbot,1); 

Dy(kscatbot+1,1)=Dy(kscatbot+1,1)-D_E(kscatbot+1)*Hxinc(kscatbot,1); 

Dx(kscattop+1,1)=Dx(kscattop+1,1)-D_E(kscattop+1)*Hyinc(kscattop+1,1); 
Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1) ; 
Ex(ze,1)=epsinvxx(ze).*Dx(ze,1)+epsinvxy(ze).*Dy(ze,1)+epsinvxz(ze).*Dz 
(ze,1); 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze, 1)+epsinvyy (ze) .*Dy(ze, 1)+epsinvyz (ze) . *Dz 
(ze,1) ; 

Ez(ze,1)=epsinvxz(ze) .*Dx(ze,1)+epsinvyz(ze) .*Dy(ze,1)+epsinvzz (ze) .*Dz 
(ze,1) ; 

Exinc(ze,1)=C_E (ze) .*Exinc(ze,1)-D_E (ze) . * (Hyinc(ze,1)-Hyinc(ze-1,1)); 
Eyinc(ze,1)=C_E(ze) .*Eyinc(ze, 1)+D_E(ze) .*(Hxinc(ze,1)-Hxinc(ze-1,1)); 
Exinc(kincstart+1,1)=-E0*sin(omega*Tbreak*dt)*(1-exp(-(Tbreak*dt- 
tdelay) A 2/sigmadelay)) ; 

Hx ( : , 2) =C_H. *Hx ( : , 1) +D_H. * (Ey (zh, 1) -Ey (zh-1, 1) ) ; 

Hy(:,2)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1,1)); 

Hx(kscatbot,2)=Hx(kscatbot,2)-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,2)=Hy(kscatbot,2)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,2)=Hx(kscattop+1,2)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 
Hy(kscattop+1,2)=Hy(kscattop+1,2)-D_H(kscattop+1)*Exinc(kscattop+1, 1) ; 
Hxinc(:,2)=C_H.*Hxinc(:,1)+D_H.* (Eyinc(zh,1)-Eyinc(zh-1,1)); 

Hyinc(:,2)=C_H.*Hyinc(:,1)-D_H.* (Exinc(zh,1)-Exinc(zh-1,1)); 

Hxinctilde(:,1)=(Hxinc(ze,1)+Hxinc(ze,2)+Hxinc(ze-1,1)+Hxinc(ze- 
1,2))/4; 

Hyinctilde(: , 1) = (Hyinc(ze,1)+Hyinc(ze,2)+Hyinc(ze-1,1)+Hyinc (ze- 
1,2))/4; 

Hxtilde(: , 1) = (Hx(ze,1)+Hx(ze, 2 ) +Hx(ze-1,1)+Hx(ze-l,2))/4; 

Hytilde(:,l) = (Hy(ze,l)+Hy(ze,2)+Hy(ze-1,1)+Hy (ze-1,2))/4; 

Poyntingx(: , 1)=-Ez(ze, 1) .*Hytilde(:,1); 

Poyntingy(:,1)=Ez(ze, 1) .*Hxtilde(:,1); 

Poyntingz(:,1)=Ex(ze,1).*Hytilde(:,1)-Ey(ze,1).*Hxtilde(:,1); 
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PoyntingMag(:,1)=sqrt(Poyntingx(: , 1) .^2+Poyntingy(: , 1) . A 2+Poyntingz ( :,1 
). A 2); 

%for n=Tbreak:Tfin 

%makes Dx(:,l) actually at n=Tfin-l-2*timeperiod 
for n=l:2*timeperiod 

Dx(ze,n+1)=C_E(ze).*Dx(ze,n)-D_E(ze). *(Hy(ze,n+1)-Hy(ze-1,n+1)); 

Dy(ze,n+1)=C_E(ze).*Dy(ze,n)+D_E(ze).*(Hx(ze,n+1)-Hx(ze-1,n+1)); 

Dx(kscatbot + 1, n+1)=Dx(kscatbot + 1,n+1)+D_E(kscatbot + 1)*Hyinc(kscatbot,n+ 

i); 

Dy (kscatbot + 1,n+1)=Dy(kscatbot + 1, n+1)- 
D_E(kscatbot + 1)*Hxinc(kscatbot, n+1); 

Dx(kscattop+1,n+1)=Dx(kscattop+1, n+1) - 
D_E(kscattop+1)*Hyinc(kscattop+1, n+1) ; 

Dy(kscattop+1,n+1)=Dy(kscattop+1,n+1)+D_E(kscattop+1)*Hxinc(kscattop+1, 
n+1) ; 

Ex(ze, n+1)=epsinvxx(ze) .*Dx(ze,n+1)+epsinvxy(ze) .*Dy(ze,n+1)+epsinvxz(z 
e) . *Dz (ze,n+1); 

Ey(ze,n+1)=epsinvxy(ze) .*Dx(ze,n+1)+epsinvyy(ze) .*Dy(ze, n+1)+epsinvyz(z 
e).*Dz(ze,n+1); 

Ez(ze, n+1)=epsinvxz(ze) .*Dx(ze, n+1)+epsinvyz(ze) .*Dy(ze,n+1)+epsinvzz (z 
e).*Dz(ze,n+1); 

Exinc(ze,n+1)=C_E (ze) .*Exinc(ze,n)-D_E(ze) . *(Hyinc(ze, n+1)-Hyinc(ze- 
1,n+1)); 

Eyinc(ze,n+1)=C_E (ze) .*Eyinc(ze,n)+D_E (ze) . *(Hxinc(ze, n+1)-Hxinc (ze- 
1,n+1)); 

Exinc(kincstart + 1, n+1)=-E0*sin(omega*(Tbreak+n)*dt)*(1-exp(- 
( (Tbreak+n)*dt-tdelay) A 2/sigmadelay)); 

Hx(:,n+2)=C_H.*Hx(:,n+1)+D_H.*(Ey(zh, n+1)-Ey(zh-1, n+1)) ; 

Hy(:,n+2)=C_H.*Hy(:,n+1)-D_H.*(Ex(zh,n+1)-Ex(zh-1, n+1)) ; 

Hx(kscatbot,n+2)=Hx(kscatbot, n+2)- 
D_H(kscatbot)*Eyinc(kscatbot+1,n+1); 

Hy(kscatbot,n+2)=Hy(kscatbot,n+2)+D_H(kscatbot)*Exinc(kscatbot + 1,n+1) ; 

Hx(kscattop+1,n+2)=Hx(kscattop+1,n+2)+D_H(kscattop+1)*Eyinc(kscattop+1, 
n+1) ; 

Hy(kscattop+1,n+2)=Hy(kscattop+1, n+2)- 
D_H(kscattop+1)*Exinc(kscattop+1, n+1) ; 

Hxinc(:,n+2)=C_H.*Hxinc(:,n+1)+D_H.*(Eyinc(zh,n+1)-Eyinc(zh-1, n+1)) ; 
Hyinc(:,n+2)=C_H.*Hyinc(:,n+1)-D_H.*(Exinc(zh,n+1)-Exinc(zh-1, n+1)) ; 
Hxinctilde(: , n+1) = (Hxinc(ze,n+1)+Hxinc(ze,n+2)+Hxinc (ze- 
1,n+1)+Hxinc(ze-1,n+2))/4; 

Hyinctilde(:,n+1)=(Hyinc(ze,n+1)+Hyinc(ze,n+2)+Hyinc(ze- 
1,n+1)+Hyinc(ze-1,n+2))/4; 

Hxtilde(:,n+1) = (Hx(ze,n+1)+Hx(ze,n+2)+Hx(ze-1,n+1)+Hx(ze-1, n+2 ))/4 ; 
Hytilde(: , n+1) = (Hy(ze,n+1)+Hy(ze,n+2)+Hy(ze-1,n+1)+Hy(ze-l,n+2))/4; 
Poyntingx(:,n+1)=-Ez(ze,n+1).*Hytilde(:,n+1); 

Poyntingy(:,n+1)=Ez(ze,n+1).*Hxtilde(:,n+1); 

Poyntingz(:,n+1)=Ex(ze,n+1).*Hytilde(:,n+1)- 
Ey(ze,n+1).*Hxtilde(:,n+1); 


93 


PoyntingMag(:,n+1)=sqrt(Poyntingx(:,n+1). A 2+Poyntingy(:,n+1). A 2+Poyntin 

gz(:,n+1) . A 2) ; 

end 

Intensity=PoyntingMag(: , 1) ; 
for n=2:2*timeperiod 

Intensity=Intensity+2*PoyntingMag(: , n) ; 

end 

Intensity=Intensity+PoyntingMag(:,2*timeperiod+l); 
Intensity=Intensity*dt/EO A 2/(4*pi/omega); 

kPMLbot5=kPMLbot; kscatbot5=kscatbot; kairbot5=kairbot; 
kglassbot5=kglassbot; 

kglasstop5=kglasstop; kairtop5=kairtop; kscattop5=kscattop; 
kPMLtop5=kPMLtop; 

Ntot5=Ntot; dz5=dz; dt5=dt; 

Ex5=Ex; Ey5=Ey; Ez5=Ez; 

Hx5=Hx; Hy5=Hy; 

Exinc5=Exinc; Eyinc5=Eyinc; 

Hxinc5=Hxinc; Hyinc5=Hyinc; 

Hxinctilde5=Hxinctilde; Hyinctilde5=Hyinctilde; 

Hxtilde5=Hxtilde; Hytilde5=Hytilde; 

Poyntingx5=Poyntingx; Poyntingy5=Poyntingy; Poyntingz5=Poyntingz; 

PoynitngMag5=PoyntingMag; Intensity5=Intensity; 

timeperiod5=timeperiod; Tfin5=Tfin; Tbreak5=Tbreak; 

kincstart5=kincstart; 

zTotE5=dz*((1:Ntot+l)-kglassbot-1); 

zTotH5=dz*((1:Ntot)-1/2-kglassbot); 

tic 


toe 


%% Finer grid 
refinement=2; 

NLCP=NLCP*refinement %number of numerical cells in the LCP 

dz=h/NLCP; 


Nglass=Nglass*refinement; 

%number 

of 

numerical 

cells 

in 

the 

glass 

Nair=Nair*refinement; 

%number 

of 

numerical 

cells 

in 

the 

air 

between glass and scattering 
Nscat=Nscat*refinement; 

layer 

%number 

of 

numerical 

cells 

in 

the 


scattering later 
NPML=NPML*refinement; 

%number 

of 

numerical 

cells 

in 

the 



perfectly matched layer 
cdtoverdz=l; 
dt=dz*cdtoverdz; 

Tfin=floor(EndTime/dt)+1; 
timeperiod=round(2*pi/omega/dt); 

sigmamax=2/dt; %10 A 17; %in the PML, sigmaz=sigmamax((kPML- 

1)/(NPML-1)) A sigmam 
sigmam=3; 

kPMLbot=NPML; % sigma defined but equal to zero at 

k=kPMLbot 


kscatbot=kPMLbot+Nscat; % k=kscatbot is considered to be in the 

total field; kscatbot-1/2 is in the scattered field 
kairbot=kscatbot+Nair; % k=kairbot is the bottom edge of the 

bottom plate 
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kglassbot=kairbot+Nglass; 
plate/where the anchoring 
kglasstop=kglassbot+NLCP; 
plate/where the anchoring 
kairtop=kglasstop+Nglass; 
plate 

kscattop=kairtop+Nair; 
kscattop+1/2 is not 
kPMLtop=kscattop+Nscat; 
k=kPMLtop 
Ntot=kPMLtop+NPML; 
ze=2:Ntot; 

in the D and E equations 
zh=2:Ntot + 1; 

z6=linspace(0,h,NLCP+1); 
kincstart=refinement*(kinc 


% k=kglassbot is the top edge of the bottom 
condition is applied 

% k=kglasstop is the bottom edge of the top 
condition is applied 

% k=kairtop is the top edge of the top 

% k=kscattop is in the total field; 

% sigma defined but equal to zero at 

% Ntot cells in total% ze=2:Ntot; 

% k values for updating the gradients of H 


%z values used for computing the orientation 
start5-kglassbot5)tkglassbot; 


%% compute the permittivity 
solinit=bvpinit(z6,@LEodeinit) ; 
sol=bvp4c(@LEodefun,@LEbcfun, solinit) ; 
u=deval(sol,z6) ; 
psiLE=u (1, : ); 
thetaLE=u (3, : ); 


nl=cos(thetaLE).*cos(psiLE); 
n2=cos(thetaLE).*sin(psiLE); 
n3=sin(thetaLE); 


epsxx=zeros(Ntot + 1, 1) ; 
epsyy=zeros(Ntot + 1,1) ; 
epszz=zeros(Ntot+1,1); 
epsxy=zeros(Ntot + 1, 1) ; 
epsxz=zeros(Ntot+1,1); 
epsyz=zeros(Ntot + 1, 1) ; 

epsxx(1:kPMLbot + 1)=1; epsyy(1:kPMLbot + 1)=1; epszz(1:kPMLbot+1)=1; 

epsxx(kPMLbot+2:kairbot)=epsilon_air; 

epsyy(kPMLbot+2:kairbot)=epsilon_air; 

epszz(kPMLbot+2:kairbot)=epsilon_air; 

epsxx(kairbot+1:kglassbot)=epsilon_glass; 

epsyy(kairbot + 1:kglasshot)=epsilon_glass; 

epszz(kairbot + 1:kglassbot)=epsilon_glass; 

epsxx(kglassbot + 1:kglasstop+1)=epsilon_perp+delta_epsilon*nl. A 2; 

epsyy(kglassbot + 1:kglasstop+1)=epsilon_perp+delta_epsilon*n2. A 2; 

epszz(kglassbot+1:kglasstop+1)=epsilon_perp+delta_epsilon*n3. A 2; 

epsxy(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n2; 

epsxz(kglassbot+1:kglasstop+1)=delta_epsilon*nl.*n3; 

epsyz(kglassbot+1:kglasstop+1)=delta_epsilon*n2.*n3; 

epsxx(kglasstop+2:kairtop+1)=epsilon_glass; 

epsyy(kglasstop+2:kairtop+1)=epsilon_glass; 

epszz(kglasstop+2:kairtop+1)=epsilon_glass; 

epsxx(kairtop+2:kPMLtop)=epsilon_air; 

epsyy(kairtop+2:kPMLtop)=epsilon_air; 

epszz(kairtop+2:kPMLtop)=epsilon_air; 
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epsxx(kPMLtop+1:Ntot)=1; epsyy(kPMLtop+1:Ntot)= 1; 
epszz(kPMLtop+1:Ntot)= 1; 


epsdet=epsxx.*(epsyy.*epszz-epsyz. A 2)-epsxy.*(epsxy.*epszz- 

epsxz.*epsyz)+epsxz.*(epsxy.*epsyz-epsyy.*epsxz); 

epsinvxx=(epsyy.*epszz-epsyz. A 2) ./epsdet; 

epsinvxy=(epsxz.*epsyz-epsxy.*epszz)./epsdet; 

epsinvxz=(epsxy.*epsyz-epsyy.*epsxz)./epsdet; 

epsinvyy=(epsxx.*epszz-epsxz. A 2)./epsdet; 

epsinvyz=(epsxy.*epsxz-epsxx.*epsyz)./epsdet; 

epsinvzz=(epsxx.*epsyy-epsxy. A 2)./epsdet; 


%% Define coefficients 

C_E=ones(Ntottl,1); 

update equations 

D_E=cdtoverdz*ones(Ntottl, 1) ; 

update equations 

C_H=ones(Ntot,1); 

update equations 

D_H=cdtoverdz*ones(Ntot, 1) ; 

update equations 

for k=l:kPMLbot 

PML 


% factor of the D A n term in the D 

% factor of the curl H term in the D 

% factor of the H A n term in the H 

% factor of the curl H term in the H 

% fill in the effect of sigma in the 


sigmaz=sigmamax*((kPMLbot-k+1)/kPMLbot) A sigmam; 
C_E(k) = (2-dt*sigmaz)/(2+dt*sigmaz) ; 

D_E (k)=2*cdtoverdz/(2+dt*sigmaz); 

C_E(Ntot+2-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_E(Ntot+2-k)=2*cdtoverdz/(2+dt*sigmaz); 
sigmaz=sigmamax*((kPMLbot-k+0.5)/kPMLbot) A sigmam; 
C_H(k) = (2-dt*sigmaz)/(2+dt*sigmaz) ; 

D_H(k)=2*cdtoverdz/(2+dt*sigmaz ) ; 

C_H(Ntottl-k)=(2-dt*sigmaz)/(2+dt*sigmaz); 

D_H(Ntottl-k)=2*cdtoverdz/(2+dt*sigmaz); 

end 


%% Define electromagnetic field variables 

Dx=zeros(Ntot+l,2*timeperiod+l); % D and E are known on the space 

time grid 

Dy=zeros(Ntot+l,2*timeperiod+l); % D and E fields are zero 

initially (n=0), but this is not stored. 


Dz=zeros(Ntot+l,2*timeperiod+l); 
k=Ntot+l but they are stored so that 
Ex=zeros(Ntot+l,2*timeperiod+l); 
z=(k-kglassbot—1)*dz, 

Ey=zeros(Ntot + l,2*timeperiod+l) ; 

Ex(k+1,n) 

Ez=zeros(Ntot + l,2*timeperiod+l) ; 
Hx=zeros(Ntot,2*timeperiod+2) ; 
z= (k-kglassbot + 1/2)*dz 
Hy=zeros(Ntot,2*timeperiod+2 ) ; 

Hx(k+1,n+1) 

Exinc=zeros(Ntot + l, 2*timeperiod+l) ; 
can be nonzero, dependending on epsi 
Eyinc=zeros(Ntot + l,2*timeperiod+l) ; 


% Dx and Ex are zero at k=l and 

% Dx(k,n) is Dx at t=n*dt and 

% or Ex_k A n is stored as 

% H is known on the half-grid. 

% Hx(k,n) is Hx at t=(n-1/2)*dt, 

% Hx_(k+1/2) A (n+1/2) is stored as 

% Hz is zero and not stored. Ez 
on 
"6 
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Hxinc=zeros(Ntot,2*timeperiod+2); % Exinc, Eyinc, Hxinc, and Hyinc 

are the incident field 
Hyinc=zeros(Ntot, 2*timeperiod+2) ; 

Hxinctilde=zeros(Ntot-1,2*timeperiod+l); % Hxinctilde interpolates 

Hxinc onto the grid 

Hyinctilde=zeros(Ntot-1,2*timeperiod+l); 

Hxtilde=zeros(Ntot-1,2*timeperiod+l); % Hxtilde interpolates Hx 

onto the grid 

Hytilde=zeros(Ntot-1,2*timeperiod+l); 

Poyntingx=zeros(Ntot-1,2*timeperiod+l); % The Poynting vector is known 

on the space-time grid, and so H must be 

Poyntingy=zeros(Ntot-1,2*timeperiod+l); % interpolated in time and in 

space. Thus it is not defined at k=l or k=Ntot+l 
Poyntingz=zeros(Ntot-1,2*timeperiod+l); % or at t=Tfin 

PoyntingMag=zeros(Ntot-1,2*timeperiod+l); 

Tbreak=Tfin-2*timeperiod-l; 

% The starting main loop, not storing every step 
Exinc(kincstart+1,1)=-E0*sin(omega*dt)*(1-exp(-(dt- 
tdelay)^2/sigmadelay)); 

Hx(kscatbot,1)=-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,1)=D_H(kscatbot)*Exinc(kscatbot + 1,1) ; 

Hxinc(: , 1)=D_H .* (Eyinc(zh,1)-Eyinc(zh-1,1)); 

Hyinc(:,1)=-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)); 
for n=2:Tbreak-l 

Dx (ze, 1) =C_E (ze) . *Dx (ze, 1) -D_E (ze) . * (Hy (ze, 1) -Hy (ze-1, 1) ) ; 

Dy(ze,1)=C_E(ze).*Dy(ze,1)+D_E(ze).*(Hx(ze,1)-Hx(ze-1,1)); 

Dx(kscatbot+1,1)=Dx(kscatbot+1,1)+D_E(kscatbot+1)*Hyinc(kscatbot,1); 
Dy(kscatbot + 1,1)=Dy(kscatbot + 1,1)-D_E(kscatbot + 1)*Hxinc(kscatbot, 1) ; 
Dx(kscattop+1,1)=Dx(kscattop+1, 1)- 
D_E(kscattop+1)*Hyinc(kscattop+1,1) ; 

Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1, 1) ; 

Ex(ze,1)=epsinvxx(ze) .*Dx(ze,1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz (ze) .*Dz 
(ze,1); 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze, 1)+epsinvyy (ze) .*Dy(ze, 1)+epsinvyz (ze) . *Dz 
(ze,1); 

Exinc(ze,1)=C_E (ze) .*Exinc(ze,1)-D_E(ze) .*(Hyinc(ze,1)-Hyinc (ze- 

i,D); 

Eyinc(ze,1)=C_E(ze).*Eyinc(ze,1)+D_E(ze).*(Hxinc(ze,1)-Hxinc(ze- 

i,D); 

Exinc(kincstart+1,1)=-E0*sin(omega*n*dt)*(1-exp(-(n*dt- 
tdelay)^2/sigmadelay)) ; 

Hx(:,1)=C_H.*Hx(:,1)+D_H.*(Ey(zh,1)-Ey(zh-1,1)); 

Hy(:,1)=C_H.*Hy(:,1)-D_H.*(Ex(zh, 1)-Ex(zh-1,1)) ; 

Hx(kscatbot,1)=Hx(kscatbot,1)-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,1)=Hy(kscatbot,1)+D_H(kscatbot)*Exinc(kscatbot + 1,1) ; 

Hx(kscattop+1,1)=Hx(kscattop+1,1)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 

Hy(kscattop+1,1)=Hy(kscattop+1,1)- 
D_H(kscattop+1)*Exinc(kscattop+1, 1) ; 

Hxinc(:,1)=C_H.*Hxinc(:,1)+D_H.*(Eyinc(zh, 1)-Eyinc(zh-1,1)) ; 

Hyinc(:,1)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh, 1)-Exinc(zh-1,1)) ; 
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end 

%n=Tbreak; 

Dx(ze,1)=C_E(ze) .*Dx(ze,1)-D_E(ze) .* (Hy (ze,1)-Hy(ze-1,1)); 

Dy(ze,1)=C_E(ze).*Dy(ze,1)+D_E(ze).*(Hx(ze,1)-Hx(ze-1,1)); 

Dx(kscatbot+1,1)=Dx(kscatbot+1,1)+D_E(kscatbot+1)*Hyinc(kscatbot,1); 

Dy(kscatbot + 1,1)=Dy(kscatbot + 1,1)-D_E(kscatbot + 1)*Hxinc(kscatbot, 1) ; 

Dx(kscattop+1,1)=Dx(kscattop+1,1)-D_E(kscattop+1)*Hyinc(kscattop+1,1); 
Dy(kscattop+1,1)=Dy(kscattop+1,1)+D_E(kscattop+1)*Hxinc(kscattop+1,1); 
Ex(ze,1)=epsinvxx(ze) .*Dx(ze,1)+epsinvxy(ze) .*Dy(ze,1)+epsinvxz (ze) .*Dz 
(ze,1); 

Ey(ze, 1)=epsinvxy(ze) .*Dx(ze, 1)+epsinvyy (ze) .*Dy(ze, 1)+epsinvyz (ze) . *Dz 
(ze,1); 

Ez(ze, 1)=epsinvxz(ze) .*Dx(ze, 1)+epsinvyz (ze) .*Dy(ze, 1)+epsinvzz (ze) . *Dz 
(ze,1) ; 

Exinc(ze,1)=C_E(ze) .*Exinc(ze,1)-D_E(ze) .*(Hyinc(ze, 1)-Hyinc(ze-1, 1)) ; 
Eyinc(ze, 1)=C_E(ze) .*Eyinc(ze,1)+D_E (ze) . *(Hxinc(ze,1)-Hxinc(ze-1,1)); 
Exinc(kincstart + 1,1)=-E0*sin(omega*Tbreak*dt)*(1-exp(-(Tbreak*dt- 
tdelay)^2/sigmadelay)); 

Hx ( : ,2)=C_H.*Hx ( : , 1) +D_H. * (Ey (zh, 1) -Ey (zh-1, 1) ) ; 

Hy(:,2)=C_H.*Hy(:,1)-D_H.*(Ex(zh,1)-Ex(zh-1,1)); 

Hx(kscatbot,2)=Hx(kscatbot,2)-D_H(kscatbot)*Eyinc(kscatbot + 1, 1) ; 

Hy(kscatbot,2)=Hy(kscatbot,2)+D_H(kscatbot)*Exinc(kscatbot+1,1); 

Hx(kscattop+1,2)=Hx(kscattop+1,2)+D_H(kscattop+1)*Eyinc(kscattop+1,1); 
Hy(kscattop+1,2)=Hy(kscattop+1,2)-D_H(kscattop+1)*Exinc(kscattop+1, 1) ; 
Hxinc(:,2)=C_H.*Hxinc(:,1)+D_H.* (Eyinc(zh,1)-Eyinc(zh-1,1)); 

Hyinc(:,2)=C_H.*Hyinc(:,1)-D_H.*(Exinc(zh,1)-Exinc(zh-1,1)); 

Hxinctilde(:,1) = (Hxinc(ze,1)+Hxinc(ze, 2 ) +Hxinc(ze-1,1)+Hxinc (ze- 
1,2))/4; 

Hyinctilde(: , 1) = (Hyinc(ze, 1)+Hyinc(ze,2)+Hyinc(ze-1,1)+Hyinc(ze- 
1,2))/4; 

Hxtilde(:,1)=(Hx(ze,1)+Hx(ze, 2 ) +Hx(ze-1,1)+Hx(ze-1,2))/4; 

Hytilde(:,1) = (Hy(ze,1)+Hy(ze, 2)+Hy(ze-1, 1)+Hy(ze-l,2))/4; 

Poyntingx(:,1)=-Ez(ze,1).*Hytilde(:,1); 

Poyntingy(:,1)=Ez(ze,1).*Hxtilde(:,1); 

Poyntingz(:,1)=Ex(ze,1).*Hytilde(:,1)-Ey(ze,1).*Hxtilde(:,1); 
PoyntingMag(:,1)=sqrt(Poyntingx(: , 1) .^2+Poyntingy(: , 1) .^2+Poyntingz(:,1 
). A 2); 

%for n=Tbreak:Tfin 

%makes Dx(:,l) actually at n=Tfin-l-2*timeperiod 
for n=l:2*timeperiod 

Dx(ze,n+1)=C_E(ze).*Dx(ze,n)-D_E(ze). *(Hy(ze,n+1)-Hy(ze-1,n+1)); 

Dy(ze,n+1)=C_E(ze).*Dy(ze,n)+D_E(ze).*(Hx(ze,n+1)-Hx(ze-1,n+1)); 

Dx(kscatbot+1,n+1)=Dx(kscatbot+1,n+1)+D_E(kscatbot+1)*Hyinc(kscatbot,n+ 

i); 

Dy (kscatbot + 1,n+1)=Dy(kscatbot + 1, n+1)- 
D_E(kscatbot + 1)*Hxinc(kscatbot, n+1); 

Dx(kscattop+1,n+1)=Dx(kscattop+1, n+1)- 
D_E(kscattop+1)*Hyinc(kscattop+1, n+1) ; 

Dy(kscattop+1,n+1)=Dy(kscattop+1,n+1)+D_E(kscattop+1)*Hxinc(kscattop+1, 
n+1) ; 

Ex(ze,n+1)=epsinvxx(ze) .*Dx(ze, n+1)+epsinvxy(ze) .*Dy(ze, n+1)+epsinvxz(z 
e).*Dz(ze,n+1); 
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Ey(ze,n+1)=epsinvxy(ze).*Dx(ze,n+1)+epsinvyy(ze).*Dy(ze,n+1)+epsinvyz(z 
e).*Dz(ze,n+1); 

Ez(ze, n+1)=epsinvxz(ze) .*Dx(ze, n+1)+epsinvyz(ze) .*Dy(ze,n+1)+epsinvzz (z 
e).*Dz(ze,n+1); 

Exinc(ze, n+1)=C_E(ze) .*Exinc(ze,n)-D_E(ze) . *(Hyinc(ze,n+1)-Hyinc (ze- 
1,n+1)); 

Eyinc(ze, n+1)=C_E(ze) .*Eyinc(ze,n)+D_E(ze) . *(Hxinc(ze,n+1)-Hxinc (ze- 
1,n+1)); 

Exinc(kincstart + 1, n+1)=-EO*sin(omega*(Tbreak+n)*dt)*(1-exp(- 
( (Tbreak+n)*dt-tdelay) A 2/sigmadelay)); 

Hx(:,n+2)=C_H.*Hx(:,n+1)+D_H.*(Ey(zh, n+1)-Ey(zh-1, n+1)) ; 

Hy(:,n+2)=C_H.*Hy(:,n+1)-D_H.*(Ex(zh,n+1)-Ex(zh-1, n+1)) ; 

Hx(kscatbot,n+2)=Hx(kscatbot, n+2)- 
D_H(kscatbot)*Eyinc(kscatbot+1,n+1); 

Hy(kscatbot,n+2)=Hy(kscatbot,n+2)+D_H(kscatbot)*Exinc(kscatbot + 1,n+1) ; 

Hx(kscattop+1,n+2)=Hx(kscattop+1,n+2)+D_H(kscattop+1)*Eyinc(kscattop+1, 
n+1) ; 

Hy (kscattop+1,n+2)=Hy(kscattop+1, n+2)- 
D_H(kscattop+1)*Exinc(kscattop+1, n+1); 

Hxinc(:,n+2)=C_H.*Hxinc(:,n+1)+D_H.*(Eyinc(zh, n+1)-Eyinc(zh-1, n+1)) ; 
Hyinc(:,n+2)=C_H.*Hyinc(:,n+1)-D_H.*(Exinc(zh,n+1)-Exinc(zh-1,n+1)); 
Hxinctilde(:,n+1)=(Hxinc(ze,n+1)+Hxinc(ze,n+2)+Hxinc(ze- 
1,n+1)+Hxinc(ze-l,n+2))/4; 

Hyinctilde(: , n+1) = (Hyinc(ze,n+1)+Hyinc(ze,n+2)+Hyinc(ze- 
1,n+1)+Hyinc(ze-l,n+2))/4; 

Hxtilde(: , n+1) = (Hx(ze, n+1)+Hx(ze, n+2)+Hx(ze-1, n+1)+Hx(ze-1, n+2))/4 ; 
Hytilde(:,n+1)=(Hy(ze,n+1)+Hy(ze,n+2)+Hy(ze-1,n+1)+Hy(ze-1,n+2))/4; 
Poyntingx(:,n+1)=-Ez(ze,n+1).*Hytilde(:,n+1); 

Poyntingy(:,n+1)=Ez(ze,n+1).*Hxtilde(:,n+1); 

Poyntingz(:,n+1)=Ex(ze,n+1).*Hytilde(:,n+1)- 
Ey(ze,n+1).*Hxtilde(:,n+1); 

PoyntingMag(: , n+1)=sqrt(Poyntingx(: , n+1) . A 2+Poyntingy(:,n+1) . A 2+Poyntin 

gz(:,n+1) . A 2); 

end 

Intensity=PoyntingMag(: , 1) ; 
for n=2:2*timeperiod 

Intensity=Intensity+2*PoyntingMag(: , n) ; 

end 

Intensity=Intensity+PoyntingMag(:,2*timeperiod+l); 
Intensity=Intensity*dt/EO A 2/(4*pi/omega); 

kPMLbot6=kPMLbot; kscatbot6=kscatbot; kairbot6=kairbot; 
kglassbot6=kglassbot; 

kglasstop6=kglasstop; kairtop6=kairtop; kscattop6=kscattop; 
kPMLtop6=kPMLtop; 

Ntot6=Ntot; dz6=dz; dt6=dt; 

Ex6=Ex; Ey6=Ey; Ez6=Ez; 

Hx6=Hx; Hy6=Hy; 

Exinc6=Exinc; Eyinc6=Eyinc; 

Hxinc6=Hxinc; Hyinc6=Hyinc; 
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Hxinctilde6=Hxinctilde; Hyinctilde6=Hyinctilde; 

Hxtilde6=Hxtilde; Hytilde6=Hytilde; 

Poyntingx6=Poyntingx; Poyntingy6=Poyntingy; Poyntingz 6=Poyntingz; 

PoynitngMag6=PoyntingMag; Intensity6=Intensity; 

timeperiod6=timeperiod; Tfin6=Tfin; Tbreak6=Tbreak; 

kincstart6=kincstart; 

zTotE6=dz*((1:Ntot+l)-kglassbot-1); 

zTotH6=dz*((1:Ntot)-1/2-kglassbot); 

toe 

%% We now will compute the true analytical solution. 

% System of Equations 
ExI=E0; 

Tplot6=2*timeperiod6tl; 

Exinctrue=EO*sin(wavenumber*(zTotE6-zTotE6(kincstart6+l))- 
omega*(Tbreak6+Tplot6-l)*dt6)*(1-exp(-((Tbreak6tTplot6-l)*dt6- 
tdelay) A 2/sigmadelay)) ; 
zstart=dz*(kincstart6-kglassbot); 

ExT=4*ExI/exp(li*wavenumber*h)*((nperp*sin(psibot) A 2) / ((1+nperp) A 2*exp( 
-li*kperp*h)- (1- 

nperp) A 2*exp(li*kperp*h))+(ntheta*(cos(psibot)) A 2)/((1+ntheta) A 2*exp(- 
li*ktheta*h)-(1-ntheta) A 2*exp(li*ktheta*h))) ; 

EyT=4*ExI/exp(li*wavenumber*h) *((- 

nperp*sin(psibot)*cos(psibot))/((1+nperp) A 2*exp(-li*kperp*h)-(1- 
nperp) A 2*exp(li*kperp*h))+(ntheta*sin(psibot)*cos(psibot)/((1+ntheta) A 2 
*exp(-li*ktheta*h)-(1-ntheta) A 2*exp(li*ktheta*h)))); 

ExR=ExI*((l-ntheta A 2)*(cos(psibot)) A 2*(1-exp(-i*2*ktheta*h))/ ( (1- 
ntheta) A 2-(1+ntheta) A 2*exp(-2*i*ktheta*h)) - (1- 
nperp A 2)*(sin(psibot)) A 2*(1-exp(2*i*kperp*h))/ ( (1- 
nperp) A 2*exp(i*2*kperp*h-(1+nperp) A 2))); 

EyR=ExI*cos(psibot)*sin(psibot)*((l-ntheta A 2)*(1-exp(- 
i*2*ktheta*h))/((1-ntheta) A 2-(1+ntheta) A 2*exp(-2*i*ktheta*h))+(1- 
nperp A 2)*(1-exp(2*i*kperp*h))/((1-nperp) A 2*exp(i*2*kperp*h)- 
(1+nperp) A 2)); 

ELplusperp=2*ExI*sin(psibot)*(1+nperp)/((1-nperp) A 2*exp(li*2*kperp*h)- 
(1+nperp) A 2); 

ELminusperp=-2*ExI*sin(psibot)*(1-nperp)/((1-nperp) A 2-(1+nperp) A 2*exp(- 
li*2*kperp*h)); 

ELplustheta=-2*ExI*cos(psibot)*(1+ntheta)/((1— 
ntheta) A 2*exp(li*2*ktheta*h)-(1+ntheta) A 2) ; 

ELminustheta=2*ExI*cos(psibot)*(1-ntheta)/((1-ntheta) A 2- 
(1+ntheta) A 2*exp(-li*2*ktheta*h)) ; 

Exyexp=exp(-li*(omega*(Tbreak6+Tplot6-l)*dt6+wavenumber*zstart+pi/2)); 
Hxyexp=exp(-li*(omega*(Tbreak6+Tplot6-l)*dt6+wavenumber*zstart+pi/2)); 

% Calculate E Components 

ExTtrue=real(Exyexp*(ExT)*exp(li*wavenumber*zTotE6(kglasstop+1:kscattop 

+ 1 ))) ; 

EyTtrue=real(Exyexp*(EyT)*exp(li*wavenumber*zTotE6(kglasstop+1:kscattop 

+ 1 ))) ; 

ExRtrue=real(Exyexp*ExR*exp(li* (- 
wavenumber*(zTotE6(kscatbot+1:kglassbot))))); 

EyRtrue=real(Exyexp*EyR*exp(li* (- 
wavenumber*(zTotE6(kscatbot + 1:kglassbot))))) ; 
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ExLtrue=real((- 

sin (psibot)*(ELplusperp*exp(li*kperp*zTotE6(kglassbot + 1:kglasstop))+ELm 
inusperp*exp(-li*kperp*zTotE6(kglassbot+1:kglasstop)))... 

+cos (psibot)* (ELplustheta*exp(li*ktheta*zTotE6(kglassbot + 1:kglasstop)) + 
ELminustheta*exp(-li*ktheta*zTotE6(kglassbot+1:kglasstop))))*Exyexp); 
EyLtrue=real((cos(psibot)*(ELplusperp*exp(li*kperp*zTotE6(kglassbot+1:k 
glasstop))+ELminusperp*exp(-li*kperp*zTotE6(kglassbot+1:kglasstop))) ... 
+ sin (psibot)*(ELplustheta*exp(i*ktheta*zTotE6(kglassbot + 1:kglasstop))+E 
Lminustheta*exp(-li*ktheta*zTotE6(kglassbot+1:kglasstop))))*Exyexp); 

% Calculate H Components 
HxTtrue=real(Hxyexp* (- 

EyT)*exp(li*wavenumber*zTotE6(kglasstop+1:kscattop+1))); 

HyTtrue=real(Hxyexp*(ExT)*exp(i*wavenumber*zTotE6(kglasstop+1:kscattop+ 
1) ) ) ; 

HxRtrue=real(Hxyexp*EyR*exp (- 

li*wavenumber*(zTotE6(kscatbot + 1:kglassbot)))) ; 

HyRtrue=real(Hxyexp*(-ExR)*exp(- 
li*wavenumber*(zTotE6(kscatbot+1:kglassbot)))); 

HxLtrue=real((- 

nperp*cos (psibot))*Hxyexp*(ELplusperp*exp(i*kperp*zTotE6(kglassbot + 1:kg 
lasstop))-ELminusperp*exp(-li*kperp*zTotE6(kglassbot+1:kglasstop))) ... 

+ (- 

ntheta*sin(psibot))*Hxyexp*(ELplustheta*exp(li*ktheta*zTotE6(kglassbot+ 
1:kglasstop))-ELminustheta*exp(- 
li*ktheta*zTotE6(kglassbot+1:kglasstop)))); 

HyLtrue=real((- 

nperp*sin(psibot))*Hxyexp*(ELplusperp*exp(li*kperp*zTotE6(kglassbot+1:k 
glasstop))-ELminusperp*exp(-li*kperp*zTotE6(kglassbot+1:kglasstop))) ... 
+(ntheta*cos(psibot))*Hxyexp*(ELplustheta*exp(i*ktheta*zTotE6(kglassbot 
+ 1:kglasstop))-ELminustheta*exp (- 
i*ktheta*zTotE6(kglassbot + 1:kglasstop)))) ; 

% Calculate E TRUE using E components 

ExTRUE=[Exinctrue(kscatbot+1:kglassbot)+ExRtrue ExLtrue ExTtrue]; 
Eyinctrue=zeros(1,length(Exinctrue)) ; 

EyTRUE=[Eyinctrue(kscatbot + 1:kglassbot)+EyRtrue EyLtrue EyTtrue] ; 

% Calculate H TRUE using H components 
Hyinctildetrue=EO*sin(wavenumber*((zTotH6(1:Ntot6- 

1)+zTotH6(2:Ntot6))/2-zTotE6(kincstart6+l))-omega*(Tbreak6+Tplot6- 
1)*dt6)*(1-exp(-((Tbreak6+Tplot6-l)*dt6-tdelay) A 2/sigmadelay)); 

HyTRUE=[Hyinctildetrue(kscatbot+1:kglassbot)+HyRtrue HyLtrue HyTtrue]; 
Hxinctildetrue=zeros(1,length(Hyinctildetrue)); 

HxTRUE=[Hxinctildetrue(kscatbot + 1:kglassbot)+HxRtrue HxLtrue HxTtrue] ; 

% Calculate Intensity TRUE 
IntensityTTRUE=zeros(size(ExTRUE)) ; 

IntensityTTRUE(kglasstop-kscatbot+1:kscattop- 
kscatbot)=4*nperp^2*sin(psibot) A 2/(4*nperp*2+ (1- 

nperp^2)*2*sin(kperp*h)^2)+4*ntheta^2*cos(psibot) A 2/(4*ntheta A 2+ (1- 
ntheta A 2) A 2*sin(ktheta*h) A 2); 

%% Plots 
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% We want to plot the last time step as a function of z for each grid. 
Tplotl=2*timeperiodl+l; 

Tplot2=2*timeperiod2+l; 

Tplot3=2*timeperiod3+l; 

Tplot4=2*timeperiod4+l; 

Tplot5=2*timeperiod5+l; 

Tplot6=2*timeperiod6+l; 

% Figure 1 compares the Ex and Hy values of the incident field. The 
green 

% vertical lines mark the top and bottom of the LCP layer. The black 
% vertical lines mark the boundaries with the scattered field, 
figure (1) 
subplot (2,1,1) 

plot(zTotEl,Exincl(:,Tplotl), 'b' ) 
hold on 

plot(zTotE2,Exinc2(:,Tplot2),' g' ) 
plot(zTotE3,Exinc3(:,Tplot3) , 'r') 
plot (zTotE4,Exinc4(:,Tplot4) , ' c ' ) 
plot(zTotE5,Exinc5(:,Tplot5) , 'k' ) 
plot (zTotE6,Exinc6(:,Tplot6), ' b' ) 
plot (zTotE6,Exinctrue, 'm: ' ) 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot],[-1.5 1.5], 'g') 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot] , [-1.5 1.5], ’ g' ) 
plot(dz*[kscatbot-kglassbot kscatbot-kglassbot], [-1.5 1.5], 'k') 
plot(dz*[kscattop-kglassbot kscattop-kglassbot], [-1.5 1.5], 'k') 
hold off 

xlim([zTotEl(1) zTotEl(Ntotl+1)]) 
title( 'Ex inc' ) 

subplot(2,1,2) 

plot((zTotHI(1:Ntotl-l)+zTotHl(2:Ntotl))/2,Hyinctildel(:,Tplotl), 'b' ) 
hold on 

plot((zTotH2(1:Ntot2-l)+zTotH2(2:Ntot2)) / 2 ,Hyinctilde2(:,Tplot2),' g' ) 
plot((zTotH3(1:Ntot3-l)tzTotH3(2:Ntot3))/2,Hyinctilde3(: ,Tplot3) , 'r ' ) 
plot((zTotH4(1:Ntot4-l)+zTotH4(2:Ntot4))/2,Hyinctilde4(:,Tplot4) , ' c' ) 
plot((zTotH5(1:Ntot5-l)tzTotH5(2:Ntot5))/2,Hyinctilde5(:,Tplot5), 'k' ) 
plot ( (zTotH6(1:Ntot6-l)tzTotH6(2:Ntot6))/2,Hyinctilde6(:,Tplot6), ' b' ) 
plot((zTotH6(l:Ntot6-l)+zTotH6(2:Ntot6))/2,Hyinctildetrue, 'm' ) 
plot(dz*[kglassbot-kglassbot kglassbot-kglassbot],[-1.5 1.5], 'g') 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot], [-1.5 1.5], 'g') 
plot(dz*[kscatbot-kglassbot kscatbot-kglassbot], [-1.5 1.5], 'k') 
plot(dz*[kscattop-kglassbot kscattop-kglassbot] , [-1.5 1.5], 'k') 
hold off 
xlabel ( ' z ' ) 

xlim([zTotHl (1) zTotHl(Ntotl-1)]) 
title ('Hy inc (averaged) ') 

% Figure 2 shows the point-wise error compared to the true solution for 
% each grid 
figure (2) 
subplot (2,1,1) 

semilogy(zTotEl,abs(Exinctrue(1:32:Ntot6+1)'-Exincl(:,Tplotl)),' b' ) 
hold on 

semilogy(zTotE2,abs(Exinctrue(1:16:Ntot6+1)'-Exinc2(:,Tplot2)),’ g' ) 
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semilogy(zTotE3, abs(Exinctrue(1:8:Ntot6 + l) '-Exinc3(: , Tplot3)) , 'r' ) 
semilogy(zTotE4,abs(Exinctrue(1:4:Ntot6+1)'-Exinc4(:,Tplot4)), 'c' ) 
semilogy(zTotE5,abs(Exinctrue(1:2:Ntot6+1)'-Exinc5(:,Tplot5)), 'k' ) 
semilogy(zTotE6,abs(Exinctrue'-Exinc6(:,Tplot6)) , ' b' ) 
hold off 

xlim([zTotEl(kscatbotl+1) zTotEl(kscattopl+1)]) 
title ('Ex inc error to true solution') 
ylim([10 A -10 10 A 0]) 

subplot(2,1,2) 

semilogy(zTotEl(2:Ntotl),abs(E0*sin(wavenumber*(zTotEl(2:Ntotl)- 
zTotEl(kincstart1+1))'-omega*(Tbreakl+Tplotl-1)*dtl)*(1-exp(- 
( (Tbreakl+Tplotl-1)*dtl-tdelay) A 2/sigmadelay))- 
Hyinctildel(:,Tplot1)) , 'b' ) 
hold on 

semilogy(zTotE2(2:Ntot2),abs(E0*sin(wavenumber*(zTotE2(2:Ntot2)- 
zTotE2(kincstart2+1))'-omega*(Tbreak2+Tplot2-l)*dt2)*(1-exp(- 
( (Tbreak2+Tplot2-l)*dt2-tdelay) A 2/sigmadelay))- 
Hyinctilde2(:,Tplot2)) , 'g' ) 

semilogy(zTotE3(2:Ntot3),abs(E0*sin(wavenumber*(zTotE3(2:Ntot3)- 
zTotE3 (kincstart3 + l)) '-omega*(Tbreak3+Tplot3-l)*dt3)*(1-exp(- 
( (Tbreak3+Tplot3-l)*dt3-tdelay) A 2/sigmadelay))- 
Hyinctilde3(:,Tplot3) ) , ' r ' ) 

semilogy(zTotE4(2:Ntot4),abs(E0*sin(wavenumber*(zTotE4(2:Ntot4)- 
zTotE4(kincstart4+1))'-omega*(Tbreak4+Tplot4-l)*dt4)*(1-exp(- 
( (Tbreak4+Tplot4-l)*dt4-tdelay) A 2/sigmadelay))- 
Hyinctilde4(: , Tplot4)) , ' c ' ) 

semilogy(zTotE5(2:Ntot5),abs(E0*sin(wavenumber*(zTotE5(2:Ntot5)- 
zTotE5(kincstart5+l))'-omega*(Tbreak5+Tplot5-l)*dt5)*(1-exp(- 
( (Tbreak5+Tplot5-l)*dt5-tdelay) A 2/sigmadelay))- 
Hyinctilde5(:,Tplot5)) , 'k' ) 

semilogy(zTotE6(2:Ntot6),abs(E0*sin(wavenumber*(zTotE6(2:Ntot6)- 
zTotE6(kincstart6+1))'-omega*(Tbreak6+Tplot6-1)*dt6)*(1-exp(- 
( (Tbreak6+Tplot6-l)*dt6-tdelay) A 2/sigmadelay))- 
Hyinctilde6(:,Tplot6) ) , 'b ' ) 
hold off 

xlim([zTotEl(kscatbotl + 1) zTotEl (kscattopl + 1)]) 
title( 'Hy inc error to true solution') 
ylim([10 A -10 10 A 0]) 
xlabel( 'z' ) 

% Figure 3 (not used) shows the the grids compared to the finest grid 
solution 

% Figure 4 shows the solutions for the full system on the different 

grids 

figure(4) 

subplot(4,1,1) 

plot (zTotEl,Exl(: , Tplotl) , ' b ' ) 
hold on 

plot(zTotE2,Ex2(:,Tplot2), 'g' ) 
plot(zTotE3,Ex3(:,Tplot3), 'r' ) 
plot (zTotE4,Ex4(:,Tplot4) , 'c' ) 
plot (zTotE5,Ex5(:,Tplot5) , 'k ' ) 
plot(zTotE6,Ex6(:,Tplot6) , 'b' ) 
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plot(zTotE6(kscatbot+1:kscattop+1),ExTRUE,' m:' ) 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot] , [-1.5 1.5], ' g' ) 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot], [-1.5 1.5], ' g' ) 
plot(dz*[kscatbot-kglassbot kscatbot-kglassbot], [-1.5 1.5], 'k') 
plot(dz*[kscattop-kglassbot kscattop-kglassbot] , [-1.5 1.5], 'k') 
hold off 

xlim([zTotEl(1) zTotEl(Ntotl+1)]) 
ylim([-1.5 1.5]) 
title (' Ex' ) 

subplot(4,1,2) 

plot (zTotEl,Eyl(:,Tplotl) , 'b' ) 
hold on 

plot(zTotE2,Ey2(:,Tplot2), 'g') 
plot(zTotE3,Ey3(:,Tplot3), 'r' ) 
plot (zTotE4,Ey4(:,Tplot4) , 'c' ) 
plot (zTotE5,Ey5(:,Tplot5) , 'k' ) 
plot(zTotE6,Ey6(:,Tplot6), 'b' ) 

plot (zTotE6(kscatbot + 1:kscattop+1),EyTRUE, ' m: ' ) 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot],[-1.5 1.5], 'g') 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot] ,[-1.5 1.5], 'g') 
plot(dz*[kscatbot-kglassbot kscatbot-kglassbot], [-1.5 1.5], 'k') 
plot(dz*[kscattop-kglassbot kscattop-kglassbot], [-1.5 1.5], 'k') 
hold off 

xlim([zTotEl(1) zTotEl(Ntotl+1)]) 
ylim([-1.5 1.5]) 
title ('Ey') 

subplot(4,1,3) % H components plotted on the x - grid 

plot((zTotHl(1:Ntotl-l)+zTotHl(2:Ntotl))/2,Hxtildel(1:Ntotl- 
1,Tplotl), 'b' ) 
hold on 

plot((zTotH2(1:Ntot2-l)+zTotH2(2:Ntot2))/2,Hxtilde2(1:Ntot2- 
1,Tplot2), 'g') 

plot((zTotH3(l:Ntot3-l)+zTotH3(2:Ntot3))/2,Hxtilde3(1:Ntot3- 
1,Tplot3), 'r' ) 

plot((zTotH4(1:Ntot4-l)+zTotH4(2:Ntot4))/2,Hxtilde4(1:Ntot4- 
1,Tplot4), 'c' ) 

plot((zTotH5(l:Ntot5-l)+zTotH5(2:Ntot5))/2,Hxtilde5(1:Ntot5- 
1,Tplot5), 'k' ) 

plot((zTotH6(l:Ntot6-l)+zTotH6(2:Ntot6))/2,Hxtilde6(1:Ntot6- 
1,Tplot6), 'b' ) 

plot(zTotE6(kscatbot+1:kscattop+1),HxTRUE,' m:' ) 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot] ,[-1.5 1.5], 'g') 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot],[-1.5 1.5], 'g') 
plot(dz*[kscatbot-kglassbot kscatbot-kglassbot], [-1.5 1.5], 'k') 
plot(dz*[kscattop-kglassbot kscattop-kglassbot], [-1.5 1.5], 'k') 
hold off 
ylim([-1.5 1.5]) 

xlim([zTotHl (1) zTotHl(Ntotl-1)]) 
title( 'Hx' ) 

subplot (4,1,4) 

plot((zTotHl(1:Ntotl-1)+zTotHl(2:Ntotl))/2,Hytildel(1:Ntotl- 
1,Tplotl), 'b' ) 
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hold on 

plot((zTotH2(1:Ntot2-l)+zTotH2(2:Ntot2))/2,Hytilde2(1:Ntot2- 
1,Tplot2), 'g' ) 

plot((zTotH3(l:Ntot3-l)+zTotH3(2:Ntot3))/2,Hytilde3(1:Ntot3- 
1,Tplot3), 'r' ) 

plot((zTotH4(1:Ntot4-l)+zTotH4(2:Ntot4))/2,Hytilde4(1:Ntot4- 
1,Tplot4), 'c' ) 

plot((zTotH5(l:Ntot5-l)+zTotH5(2:Ntot5))/2,Hytilde5(1:Ntot5- 
1,Tplot5), 'k' ) 

plot((zTotH6(l:Ntot6-l)+zTotH6(2:Ntot6))/2,Hytilde6(1:Ntot6- 
1,Tplot6), 'b') 

plot(zTotE6(kscatbot+1:kscattop+1),HyTRUE,' m:' ) 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot] , [-1.5 1.5], ’ g' ) 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot],[-1.5 1.5], 'g') 
plot(dz*[kscatbot-kglassbot kscatbot-kglassbot], [-1.5 1.5], 'k') 
plot(dz*[kscattop-kglassbot kscattop-kglassbot], [-1.5 1.5], 'k') 
hold off 
xlabel ( ' z ' ) 
title (' Hy' ) 
ylim([-1.5 1.5]) 

xlim([zTotHl (1) zTotHl(Ntotl-1)]) 

% Figure 5 shows the the grids compared to the finest grid solution 
figure (5) 
subplot (2,1,1) 

semilogy(zTotEl(kscatbotl+1:kscattopl+1),abs(ExTRUE(1:32:(kscattop- 
kscatbot + 1)) '-Exl(kscatbotl + 1:kscattopl + 1,Tplotl)) , 'b' ) 
hold on 

semilogy(zTotE2(kscatbot2 + l:kscattop2 + l) , abs(ExTRUE(1:16: (kscattop 
kscatbot + 1)) '-Ex2(kscatbot2 + l:kscattop2 + l, Tplot2)) , 'g ' ) 

semilogy(zTotE3(kscatbot3+l:kscattop3+l),abs(ExTRUE(1:8:(kscattop- 
kscatbot + 1)) '-Ex3(kscatbot3 + l:kscattop3 + l, Tplot3)) , 'r' ) 

semilogy(zTotE4(kscatbot4+l:kscattop4+l),abs(ExTRUE(1:4:(kscattop- 
kscatbot + 1)) '-Ex4(kscatbot4 + l:kscattop4 + l, Tplot4)) , 'c ' ) 

semilogy(zTotE5(kscatbot5+l:kscattop5+l),abs(ExTRUE(1:2:(kscattop- 
kscatbot+1))'-Ex5(kscatbot5+l:kscattop5+l,Tplot5)), 'k' ) 
semilogy(zTotE6(kscatbot + 1:kscattop+1) , abs(ExTRUE'- 
Ex6(kscatbot6 + l:kscattop6 + 1,Tplot6)) , 'm:' ) 
hold off 

xlim([zTotEl(kscatbotl + 1) zTotEl (kscattopl + 1)]) 
title ('Ex tilde error to Ex True') 
ylim([10 A -6 10 A 0]) 

subplot(2,1,2) 

semilogy(zTotEl(kscatbotl+1:kscattopl+1),abs(HyTRUE(1:32:(kscattop- 
kscatbot + 1)) '-Hytildel(kscatbotl + 1:kscattopl + 1, Tplotl)) , 'b' ) 
hold on 

semilogy(zTotE2(kscatbot2 + l:kscattop2 + l) , abs(HyTRUE(1:16: (kscattop 
kscatbot+1))'-Hytilde2(kscatbot2+l:kscattop2+l,Tplot2)), 'g' ) 

semilogy(zTotE3(kscatbot3+l:kscattop3+l),abs(HyTRUE(1:8:(kscattop- 
kscatbot + 1)) '-Hytilde3(kscatbot3 + l:kscattop3 + l, Tplot3)) , 'r' ) 

semilogy(zTotE4(kscatbot4+l:kscattop4+l),abs(HyTRUE(1:4:(kscattop- 
kscatbot + 1)) '-Hytilde4(kscatbot4 + l:kscattop4 + l,Tplot4)) , 'c' ) 

semilogy(zTotE5(kscatbot5+l:kscattop5+l),abs(HyTRUE(1:2:(kscattop- 
kscatbot+1))'-Hytilde5(kscatbot5+l:kscattop5+l,Tplot5)), 'k' ) 
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semilogy(zTotE6(kscatbot + 1:kscattop+1) , abs(HyTRUE'- 
Hytilde6(kscatbot + 1:kscattop+1,Tplot6)) , ' m: ' ) 
hold off 

xlim([zTotEl(kscatbotl+1) zTotEl(kscattopl+1)]) 
title( 'Hy tilde error to Hy True ') 
ylim([10 A -6 10 A 0]) 
xlabel ( ' z ' ) 

% Figure 6 shows the components of Poynting vector and the intentsity 
figure (6) 
subplot(2,2,1) 

plot((zTotHl(1:Ntotl-l)+zTotHl(2:Ntotl))/2,Poyntingxl(:,Tplotl), 'b' ) 
hold on 

plot((zTotH2(1:Ntot2-l)+zTotH2(2:Ntot2))/2,Poyntingx2(:,Tplot2), 'g' ) 
plot((zTotH3(l:Ntot3-l)+zTotH3(2:Ntot3))/2,Poyntingx3(:,Tplot3) , 'r' ) 
plot ( (zTotH4(1:Ntot4-l)+zTotH4(2:Ntot4))/2,Poyntingx4(:,Tplot4) , 'c' ) 
plot((zTotH5(l:Ntot5-l)+zTotH5(2:Ntot5))/2,Poyntingx5(:,Tplot5), 'k' ) 
plot((zTotH6(l:Ntot6-l)+zTotH6(2:Ntot6))/2,Poyntingx6(:,Tplot6), 'b' ) 
plot((zTotH6(1:Ntot6- 

l)+zTotH6(2:Ntot6))/2,PoyntingxTRUE(:,Tplot6), 'm:' ) 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot] ,[-1.25 1.25], 'g') 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot],[-1.25 1.25], 'g') 
hold off 

title (' Poyntingx ') 
ylim([-1.25 1.25]) 

xlim([zTotEl (1) zTotEl(Ntotl-1)]) 
subplot (2,2,2) 

plot((zTotHl(1:Ntotl-1)+zTotHl(2:Ntotl))/2,Poyntingyl(:,Tplotl), 'b' ) 
hold on 

plot((zTotH2(1:Ntot2-l)+zTotH2(2:Ntot2))/2,Poyntingy2(:,Tplot2), 'g' ) 
plot ( (zTotH3(l:Ntot3-l)+zTotH3(2:Ntot3))/2,Poyntingy3(:,Tplot3), 'r' ) 
plot((zTotH4(!:Ntot4-l)+zTotH4(2:Ntot4))/2,Poyntingy4(:,Tplot4), 'c' ) 
plot ( (zTotH5(l:Ntot5-l)+zTotH5(2:Ntot5))/2,Poyntingy5(:,Tplot5), 'k' ) 
plot((zTotH6(l:Ntot6-l)+zTotH6(2:Ntot6))/2,Poyntingy6(:,Tplot6) , 'b' ) 
plot((zTotH6(l:Ntot6- 

l)+zTotH6(2:Ntot6))/2, PoyntingyTRUE(: , Tplot6) , 'm: ' ) 
title (' Poyntingy ') 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot],[-1.25 1.25], 'g') 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot],[-1.25 1.25], 'g') 
hold off 

ylim([-1.25 1.25]) 

xlim([zTotEl(1) zTotEl(Ntotl-1)]) 

subplot(2,2,3) 

plot((zTotHl(1:Ntotl-1)+zTotHl(2:Ntotl))/2,Poyntingzl(:,Tplotl), 'b’ ) 
hold on 

plot((zTotH2(1:Ntot2- 

1)+zTotH2(2:Ntot2 ))/ 2 ,Poyntingz2(:,Tplot2), 'g' ) 
plot((zTotH3(1:Ntot3- 

l)+zTotH3(2:Ntot3))/2, Poyntingz3(: , Tplot3) , 'r' ) 
plot((zTotH4(1:Ntot4- 

1)+zTotH4(2:Ntot4))/2,Poyntingz4(:,Tplot4), 'c' ) 
plot((zTotH5(1:Ntot5- 

l)+zTotH5(2:Ntot5))/2, Poyntingz5(: , Tplot5) , 'k' ) 
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plot((zTotH6(1:Ntot6- 

l)+zTotH6(2:Ntot6))/2,Poyntingz6(:,Tplot6), 'b' ) 
plot((zTotH6(1:Ntot6- 

l)+zTotH6(2:Ntot6))/2, PoyntingzTRUE(: , Tplot6) , 'm: ' ) 
title (' Poyntingz ') 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot],[-1.5 1.5], 'g') 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot] , [-1.5 1.5], 'g') 
hold off 

ylim([-0.3 1.25]) 

xlim([zTotHl (1) zTotHl(Ntotl)]) 

subplot (2,2,4) 

plot((zTotHl(1:Ntotl-1)+zTotHl(2:Ntotl))/2,Intensityl,' b' ) 
hold on 

plot((zTotH2(1:Ntot2-l)+zTotH2(2:Ntot2)) 12 , Intensity2,' g' ) 
plot ( (zTotH3(l:Ntot3-l)+zTotH3(2:Ntot3))/2,Intensity3, ' r' ) 
plot ( (zTotH4(1:Ntot4-l)+zTotH4(2:Ntot4))/2,Intensity4, ' c' ) 
plot((zTotH5(1:Ntot5-l)+zTotH5(2:Ntot5))/2,Intensity5,' k' ) 
plot ( (zTotH6 (1:Ntot6-1)+zTotH6(2:Ntot6))/2,Intensity6, ’ b’ ) 
plot ( (zTotH6(l:Ntot6-l)+zTotH6(2:Ntot6))/2,IntensityTRUE, 'm: ' ) 
title( 'Intensity' ) 

plot(dz*[kglassbot-kglassbot kglassbot-kglassbot],[-1.5 1.5],' g' ) 
plot(dz*[kglasstop-kglassbot kglasstop-kglassbot],[-1.5 1.5], 'g') 
hold off 

ylim([0.38 0.55]) 

xlim([zTotEl(kscatbotl+1) zTotEl(kscattopl+1)]) 

% Figure 7 shows the error of the Poynting vector components compared 
to 

% the finest solution. The x- and y- components are identically zero 
for 

% all grids. 

figure (7) 

subplot (2,2,1) 

semilogy((zTotHl(1:Ntotl- 

1)+zTotHl(2:Ntotl))/2,abs(Poyntingx6(32:32:Ntot6-2,Tplot6)- 
Poyntingxl(:,Tplotl)) , 'b' ) 
hold on 

semilogy((zTotH2(1:Ntot2- 

1)+zTotH2(2:Ntot2))/2,abs(Poyntingx6(16:16:Ntot6-2,Tplot6)- 
Poyntingx2(:,Tplot2)) , ' g' ) 

semilogy((zTotH3(1:Ntot3- 

1)+zTotH3(2:Ntot3))/2,abs(Poyntingx6(8:8:Ntot6-2,Tplot6)- 
Poyntingx3(:,Tplot3)) , ' r' ) 

semilogy((zTotH4(1:Ntot4- 

1)+zTotH4(2:Ntot4))/2,abs(Poyntingx6(4:4:Ntot6-2, Tplot6)- 
Poyntingx4(:,Tplot4)) , 'c ' ) 

semilogy((zTotH5(1:Ntot5- 

l)+zTotH5(2:Ntot5))/2,abs(Poyntingx6(2:2:Ntot6-2,Tplot6)- 
Poyntingx5(:,Tplot5)) , 'k' ) 
hold off 

xlim([zTotEl(kairbotl + 1) zTotEl (kairtopl + 1)]) 
ylim([10 A -5 1]) 

title( ' IPoyntingx (fine) - Poyntingx (course) | ') 
subplot(2,2,2) 
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semilogy((zTotHI(1:Ntotl- 

1)+zTotHl(2:Ntotl))/2,abs(Poyntingy6(32:32:Ntot6-2,Tplot6)- 
Poyntingxl(:,Tplotl)) , 'b') 
hold on 

semilogy((zTotH2(1:Ntot2- 

1)+zTotH2(2:Ntot2))/2,abs(Poyntingy6(16:16:Ntot6-2,Tplot6)- 
Poyntingy2(:,Tplot2)) , 'g' ) 

semilogy((zTotH3(1:Ntot3- 

l)+zTotH3(2:Ntot3))/2,abs(Poyntingy6(8:8:Ntot6-2,Tplot6)- 
Poyntingy3(:,Tplot3)) , 'r' ) 

semilogy((zTotH4(1:Ntot4- 

1)+zTotH4(2:Ntot4))/2,abs(Poyntingy6(4:4:Ntot6-2,Tplot6)- 
Poyntingy4(:,Tplot4)) , 'c' ) 

semilogy((zTotH5(1:Ntot5- 

l)+zTotH5(2:Ntot5))/2,abs(Poyntingy6(2:2:Ntot6-2,Tplot6)- 
Poyntingy5(:,Tplot5)), 'k' ) 
hold off 

xlim([zTotEl(kairbotl + 1) zTotEl (kairtopl + 1)]) 

title( 'IPoyntingy (fine) - Poyntingy (course)| ') 

ylim([10 A -5 1]) 

subplot(2,2,3) 

semilogy((zTotHI(1:Ntotl- 

1)+zTotHl(2:Ntotl))/2,abs(Poyntingz6(32:32:Ntot6-2,Tplot6)- 
Poyntingzl(:,Tplotl)) , 'b’ ) 
hold on 

semilogy((zTotH2(1:Ntot2- 

1)+zTotH2(2:Ntot2))/2,abs(Poyntingz6(16:16:Ntot6-2,Tplot6)- 
Poyntingz2(:,Tplot2)) , 'g' ) 

semilogy((zTotH3(1:Ntot3- 

l)+zTotH3(2:Ntot3))/2,abs(Poyntingz6(8:8:Ntot6-2,Tplot6)- 
Poyntingz3(:,Tplot3)) , 'r' ) 

semilogy((zTotH4(1:Ntot4- 

1)+zTotH4(2:Ntot4))/2,abs(Poyntingz6(4:4:Ntot6-2,Tplot6)- 
Poyntingz4(:,Tplot4)) , 'c ' ) 

semilogy((zTotH5(1:Ntot5- 

l)+zTotH5(2:Ntot5))/2,abs(Poyntingz6(2:2:Ntot6-2,Tplot6)- 

Poyntingz5(:,Tplot5)) , 'k’ ) 

hold off 

ylim([10 A -5 1]) 

title( 'IPoyntingz (fine) - Poyntingz (course)| ') 

subplot(2,2,4) 

semilogy((zTotHI(1:Ntotl- 

1)+zTotHl(2:Ntotl))/2,abs(Intensity6(32:32:Ntot6-2)-Intensityl) , 'b’ ) 
hold on 

semilogy((zTotH2(1:Ntot2- 

1)+zTotH2(2:Ntot2))/2,abs(Intensity6(16:16:Ntot6-2)-Intensity2), 'g' ) 
semilogy((zTotH3(1:Ntot3- 

1)+zTotH3(2:Ntot3))/2,abs(Intensity6(8:8:Ntot6-2)-Intensity3),' r' ) 
semilogy((zTotH4(1:Ntot4- 

1)+zTotH4(2:Ntot4))/2,abs(Intensity6(4:4:Ntot6-2)-Intensity4),' c' ) 
semilogy((zTotH5(1:Ntot5- 

1)+zTotH5(2:Ntot5))/2,abs(Intensity6(2:2:Ntot6-2)-Intensity5),' k' ) 
hold off 

xlim([zTotEl(kairbotl+1) zTotEl(kairtopl+1)]) 
title ( ' |Intensity (fine) - Intensity (course) | ') 
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ylim([10 A -5 1]) 


% Figure 8 is the plot of the LCP director angles psi and theta, 
figure(8) 

plot(z6,psiLE, 'b') 
hold on 

plot(z6,thetaLE, 'r' ) 
hold off 

title (' Director angles \psi (blue) and \theta (red) ') 
xlabel( 'z' ) 


B. MATLAB CODE FOR CALLED FUNCTIONS 

The following functions are called by the main MATLAB program. The three 
functions, LEbcfun.m, LEodefun.m, and LEodeinit.m, were created by Dr. Eric Choate 
and not modified for puposes of this thesis. 

function res = LEbcfun(ya,yb) 

global psibot thetabot psitop thetatop; 

res = [ ya(l)-psibot 

ya(3)-thetabot 
yb(1)-psitop 
yb(3)-thetatop ]; 

function dydx=LEodefun(x,y) 
dydx=[ y (2) 

2*tan(y(3))*y(2)*y (4) 
y (4) 

-sin(2*y(3) ) *y(2) A 2/2] ; 

function yinit = LEodeinit(x) 
global h psibot thetabot psitop thetatop; 
psi=psibot+x*(psitop-psibot)/h; 
theta=thetabot+x*(thetatop-thetabot)/h; 
yinit = [ psi 

(psitop-psibot)/h 
theta 

(thetatop-thetabot)/h]; 
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